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1  INTRODUCTION 

1.1  Background  Review 

Airborne  measurements  of  gravity  were  initiated  back  in  the  late  fifties  when  the  concept 
of  airborne  gravimetry  was  first  tested  by  Thompson  (1959).  Although  the  possibility 
of  success  seemed  rather  remote  at  that  time  considerable  work  continued  in  airborne 
gravimetry  by  Nettleton  et.  al.  (1960),  Coons  et.  al.  (1962),  Gumert  and  Cobh  (1970), 
Szabo  and  Anthony  (1971),  La  Coste  et.  al.  (1977).  Only  recently,  however,  airborne 
gravity  surveys  have  met  with  considerable  success  as  reported  by  Hammer  (1982,  1983) 
and  Brozena  ( 1984).  The  major  advantages  of  airborne  gravimetry  include  speed,  efficiency, 
uniformity  of  data  quality  and  coverage  over  otherwise  inaccessible  areas.  The  major 
problem  with  airborne  gravimetry,  however,  is  the  inability  to  separate  gravitational  forces 
from  inertial  forces  ( Meissl ,  1970). 

The  problem  of  separating  inertial  and  gravitational  forces  motivated  the  development 
of  airborne  gravity  gradiometry  where  the  inertial  and  gravitational  forces  are  separated 
by  mounting  three  gravity  gradiometer  sensors  on  an  inertially  stabilized  platform  ( Moritz , 
1967,  1971).  A  gravity  gradiometer  measures  the  six  second  order  gravity  gradients.  By 
eliminating  the  effect  of  the  earth’s  normal  gravity  field  from  these  six  gravity  gradients, 
one  obtains  the  gravity  gradients  of  the  anomalous  potential.  These  second-order  gravity 
gradients  provide  short -wavelength  (high  frequency)  information  particularly  suitable  for 
a  precise  determination  of  the  anomalous  gravity  field  in  a  local  area.  In  addition,  the 
high  frequencies  of  the  anomalous  gravity  field  provide  useful  information  for  geophysical 
prospecting  ( Jordan .  197S;  Brown,  1981). 

Hardware  development  for  the  gravity  gradiometers  began  towards  the  end  of  the  six¬ 
ties  and  has  continued  since  then.  During  the  seventies  four  gravity  gradiometers  were 
under  development:  the  Rotating  Gravity  Gradiometer  of  Hughes  Research  Laboratories 
( Forward ,  1971 ),  the  Spherical  Floated  Gravity  Gradiometer  of  Charles  Stark  Draper  Lab¬ 
oratory  ( Trageser,  1970,  1975),  the  Rotating  Accelerometer  Gravity  Gradiometer  of  Bell 
Aerospace/Textron  ( Metzger  and  Jircitano,  1977,  1981)  and  the  Superconducting  Grav¬ 
ity  Gradiometer  of  the  University  of  Maryland  ( Paik ,  1976,  1981,  1985).  Of  these  four, 
only  two  are  currently  under  development,  the  Bell  Rotating  Accelerometer  Gravity  Gra- 
diometer  system  now  known  as  the  Gravity  Gradiometer  Survey  System  (GGSS)  and  the 
Superconducting  Gravity  Gradiometer  developed  by  Paik  and  his  coworkers  at  the  Uni¬ 
versity  of  Maryland.  The  Superconducting  Gravity  Gradiometer  is  planned  to  be  used  for 
satellite  gradiometry.  ( Paik  et.  al.,  197S)  The  GGSS  is  the  one  which  will  be  employed 
for  airborne  gradiometery. 

1.2  Purpose  and  Scope 

The  purpose  of  this  research  is  to  develop  a  methodology  for  processing  GGSS  survey 
data.  The  GGSS  consists  of  three  gravity  gradiometer  instruments  (GGI’s)  mounted  in  an 
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PLANE  AZIMUTH  AXIS 


Figure  1:  GGI  Orientation 


‘umbrella’  configuration  of  35  degrees  as  shown  in  Figure  1  aboard  a  stabilized  platform 
together  with  the  associated  electronics,  GPS  receiver,  and  a  Rolm  computer.  The  platform 
contains  gyros  and  accelerometers  and  operates  as  a  navigator  using  GPS  position  updates, 
thus  establishing  a  coordinate  frame  for  gradiometer  measurements.  To  define  the  complete 
gravity  gradient  field  requires  three  GGI’s.  The  umbrella  configuration  was  chosen  so  that 
each  GGI  will  be  in  the  same  inertial  environment.  A  total  of  six  gradient  measurements 
are  made  in  the  orthogonal  system  of  the  umbrella  angle,  which  is  sufficient  to  define  the 
total  gravity  field. 

Gravity  gradients  are  sensed  by  differencing  the  accelerations  measured  by  two  high 
precision  Bell  Aerospace  accelerometers.  Two  pairs  of  accelerometers  are  mounted  along 
the  edge  of  a  nine  inch  disk,  with  their  sensitive  axes  tangent.  Outputs  of  opposite  pairs, 
with  sensitive  axes  180  degrees  out  of  phase,  are  summed,  and  the  outputs  of  the  pairs 
are  differenced.  With  this  arrangement,  any  linear  accelerations  to  the  disk  will  not  be 
registered  by  the  combined  output.  However,  if  an  acceleration  field  exists  with  the  specific 
force  higher  on  one  side  of  the  disk  than  the  other  ,  an  unbalance  will  result  and  the  slope 
of  the  acceleration  field  will  be  the  output. 

In  the  Bell  Aerospace  GGI’s,  the  disk  is  slowly  rotated.  When  rotation  occurs,  the 
combined  output  becomes  sinusoidal  with  the  amplitude  proportional  to  the  magnitude 
of  the  gradient  and  the  frequency  equal  to  two  times  the  rotation  rate.  This  procedure 
has  the  effect  of  shifting  the  frequency  of  the  gradient  signal  from  DC  to  twice  the  wheel 
rate.  It  now  becomes  possible  to  bandpass  filter  the  gradient  outputs,  eliminating  much  of 
the  noise.  When  the  signal  is  demodulated,  we  get  inline  and  cross  gradients.  Thus,  the 
outputs  of  the  GGSS,  recorded  on  the  data  tape,  are  six  gravity  gradient  elements  in  the 
umbrella  angle  reference  frame. 

The  data  processing  falls  naturally  into  two  stages.  The  first  stage  is  temporal  wherein 
the  data  processing  demodulates  the  survey'  system  output,  compensates  the  output  for 
the  effects  of  pressure,  temperature,  drift  rate  and  self  gradients  and  low  pass  filters  the 
data  to  produce  a  data  tape  with  one  second  samples  of  the  gradients.  Gradients  are  still 
in  the  instrument  coordinate  system,  the  umbrella  angle.  The  second  stage  is  spatial  and, 
therefore,  two  dimensional  in  nature  and  consists  of  gridding,  terrain  correction,  interpo¬ 
lation,  smoothing,  integration,  downward  continuation,  and  incorporation  of  astrogeodetic 
tie  points.  The  compensated,  demodulated,  and  filtered  gradient  signal  after  Stage  I  pro¬ 
cessing  will  be  lined  up  with  time  synchronized  position  data  from  the  aircraft  on-board 
computer  and  sampled  spatially  at  1  Km  intervals  corresponding  to  the  survey  grid.  Ide¬ 
ally,  the  1  Km  data  obtained  from  the  Stage  I  processing  will  coincide  with  the  survey  grid. 
However,  due  to  navigation  and  aircraft  control  errors,  the  actual  tracks  flown  will  deviate 
from  the  ideal  survey  tracks  by  up  to  a  few  tenths  of  a  kilometer.  The  initial  operation 
in  Stage  II  processing  will  be  a  data  gridding  exercise  where  an  interpolation  algorithm 
fa  local  least  squares  coilocator)  will  be  used  to  line  up  the  actual  gradient  data  with  the 
survey  grid. 

Stage  I  data  processing  as  well  as  the  initial  data  gridding  of  Stage  II  axe  not  within 


3 


the  scope  of  this  effort  find  henceforth  not  discussed  in  this  report  any  further.  Hence, 
it  is  assumed  that  Stage  I  processing  and  Stage  II  data  gridding  have  been  successfully 
performed  with  the  end  product  being  regular  two-dimensional  grids  of  gravity  gradient 
measurements  of  the  anomalous  potential.  Furthermore,  it  is  assumed  that  the  altitude  of 
the  survey  airplane  is  constant  within  small  bounds. 

1.3  Problem  Statement 

The  objective  of  this  research  effort  was  to  devise  a  methodology  to  process  two-dimensional 
grids  of  gravity  gradients  at  survey  altitude  to  yield  gravity  disturbance  vector  estimates 
at  the  surface  of  the  earth. 


1.4  Research  Objective 

The  measured  gravity  gradients  are  the  six  elements  of  the  gradient  tensor.  An  actual 
airborne  gradiometer  survey  over  an  area  300x300  Km  consisting  of  bidirectional  flight 
paths  5  Km  apart  with  along-track  sampling  intervals  of  1  Km  will  result  in  approximately 
220,000  measurements.  The  main  problem  with  the  determination  of  the  gravity  field  from 
airborne  gradiometery  is  the  huge  amount  of  gradient  data  collected  during  a  gradiometry 
survey.  The  primary  objective  of  this  research  effort  was  to  solve  the  problem  of  process¬ 
ing  all  the  airborne  gravity  gradient  measurements  simultaneously  in  a  computationally 
efficient  manner  without  neglecting  gradiometer  measurement  noise. 

1.5  Survey  of  Techniques 

GGSS  data  processing  methods  proposed  to  date  are  primarily  based  on  least  squares 
collocation,  an  excellent  review  of  which  is  presented  in  Moritz  (1980).  The  method  of 
least  squares  collocation  was  employed  by  Schwarz  (1977)  for  the  computation  of  the 
vertical  gravity  disturbance  vector  from  simulated  gravity  gradient  data.  Due  to  the  limited 
extent  of  the  survey  data  a  low  frequency  part  had  been  subtracted  which  resulted  in  the 
estimation  of  the  residual  gravity  field.  His  simulation  results  showed  that  15'  x  15'  mean 
disturbance  could  be  estimated  to  an  accuracy  of  2.2  mgals  with  a  profile  spacing  of  20'. 
In  tb  "  same  study  it  was  found  that  the  estimates  wil.  get  worse  if  the  gradient  data  are 
corrupted  by  systematic  errors. 

Least  squares  collocation  involves  the  inversion  of  a  matrix  of  order  equal  to  the  number 
of  measurements.  Since  the  number  of  measurements  collected  during  the  airborne  GGSS 
test  is  expected  to  be  approximately  2  x  105,  a  full  scale  least  squares  collocation  approach 
is  presently  impractical.  A  template  method  ( Goldstein  and  White,  1985)  exploits  the 
attenuation  with  distance  in  the  correlations  between  gravity  gradients.  By  averaging  the 
measurements  over  compartments  whose  sizes  increase  with  distance  from  the  estimation 
point  the  number  and  density  of  observations  are  reduced  significantly.  Thus,  the  short 
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wavelengths  of  gradients  in  distant  zones  which  signify  local  density  contrasts  and  have 
minimal  correlation  with  the  estimate  are  averaged  out,  while  the  longer  wavelengths 
are  retained.  The  averaging  compartments  form  a  template  centered  directly  above  the 
estimation  point,  hence  the  name  'template  method’.  Once  the  data  set  with  large  numbers 
of  points  has  been  transformed  to  a  much  smaller  set  of  averages  the  gravity  disturbance 
vector  is  estimated  using  space  domain  least  squares  collocation.  This  method  requires 
generation  of  covariances  between  gravity  gradient  averages  implying  the  necessity  of  a 
covariance  model  that  is  amenable  to  analytic  integration.  The  covariance  function  used 
for  the  template  method  is  the  Attenuated  White  Noise  (AWN)  model  (Heller  and  Jordan , 
1979).  The  template  method  is  simplest  for  uniformly  spaced  data,  since  a  fixed  template 
design  can  be  utilized.  However,  for  non  uniformly  spaced  data  the  template  needs  to  be 
restructured  and  the  weights  for  each  template  zone  recomputed  again  for  each  estimation 
point.  To  redefine  the  zones  and  recompute  the  optimal  weights  a  matrix  of  order  equal 
to  the  number  of  template  zones  needs  to  be  inverted. 

Another  space-domain  method  for  processing  GGSS  data  is  based  on  Rummel’s  (19S2) 
technique  of  developing  an  integral  formula  estimator  that  is  optimized  for  the  effect  of 
measurement  noise.  This  method  was  developed  as  an  alternative  to  the  inverse  Stokes 
integral  for  processing  large,  densely  spaced  satellite  altimetry  data  without  ignoring  mea¬ 
surement  errors.  Assuming  flat  earth  approximation  Jekeli  (1984)  extended  the  idea  to 
incorporate  heterogeneous  gravity  gradient  data  sets  at.  altitude.  The  appropriate  Ker¬ 
nel  function  is  expressed  as  a  Fourier  transform  of  its  spectrum  which  in  the  flat  earth 
approximation  is  fitted  by  pieces  of  powers  of  spatial  frequency  and  thereby  amenable  to 
analytical  integration.  Thus,  from  each  set  of  measured  gradients  the  whole  gravity  distur¬ 
bance  vector  can  be  recovered.  Gravity  disturbance  vector  estimates  from  different  sets  of 
gravity  gradient  measurements  can  be  combined  in  a  least  squares  adjustment  to  provide 
the  final  estimate  ( Jekeli ,  1984).  A  further  enhancement  to  the  estimator  was  proposed  by 
Jekeli  (1985)  wherein  a  method  to  integrate  isolated  surface  tie-point,  data  was  described. 
However,  the  estimator  assumes  infinitely  extended  gradient  data  and  does  not  assign  suf¬ 
ficient  weight  to  the  tie-point  data.  Thus,  there  is  no  reproducibility  at  the  tie-point,  data 
as  with  optimal  least-squares  collocation.  Hence,  the  incorporation  of  tie-point  data  is  not 
optimal  in  this  estimator.  On  the  other  hand,  without  incorporating  tie-points  Jekeli's 
(1985)  simulation  results  show  that  a  noticeable  bias  and  trend  could  not  be  estimated 
from  the  given  gradient  data.  Also,  because  of  the  lack  of  data  outside  the  survey  area, 
estimates  near  the  edge  are  distorted  and  unreliable.  One  limitation  of  this  method  is  the 
requirement  that  the  data  be  on  a  plane  of  constant  altitude.  Furthermore,  to  recover  the 
gravity  disturbance  vector  at  the  surface  of  the  earth  downward  continuation  has  to  be 
performed. 

Another  technique  based  on  Stokes’  theorem  has  been  proposed  by  Rufiy  ( 1986)  wherein 
the  connection  between  cross  products  of  vector  fields  and  line  integrals  of  gradients  is 
exploited.  For  conservative  fields,  the  constraint  of  zero  cross  product  of  the  vector  field 
is  equivalent  to  zero  closed  line  integrals  of  the  gradients.  This  fact  is  exploited  in  this 


method  to  recognize  that  due  to  the  presence  of  noise  these  gradient  line  integrals  are 
path  dependent.  The  method  computes  the  best  estimates  of  segment  integrals  and  then 
combines  the  estimates  from  the  various  paths.  Thus,  the  algorithm  simultaneously  weights 
all  possible  path  integrals  of  the  gradient  to  yield  an  optimal  estimate  of  the  disturbance. 
The  method  considers  all  paths  simultaneously  but  does  not  include  any  modeling  of  the 
statistics  of  the  field.  The  gravity  estimates  are  obtained  at  the  survey  altitude  from 
gradient  track  data  at  altitude.  Thus,  to  recover  the  gravity  disturbance  vector  at  the 
surface  of  the  earth  downward  continuation  has  to  be  performed. 

A  combination  of  least-squares  collocation  in  the  space  domain  and  Wiener  filtering 
in  the  frequency  domain  has  been  proposed  by  Bell  Aerospace/Textron  ( Hutcheson ,  1985; 
Hutcheson  and  Grierson ,  1985).  The  least  squares  collocator  is  used  to  estimate  the  low 
frequency  component  of  gravity  disturbance  on  a  sparse  grid.  This  is  then  combined  in  the 
frequency  domain  with  the  complementary  high  frequency  disturbance  estimate  obtained 
from  the  Wiener  smoother.  The  resulting  map  is  then  transformed  back  into  the  space 
domain.  For  both  the  low  and  high  frequency  parts  of  the  spectrum,  plane  integration  at 
flying  altitude  and  downward  continuation  are  carried  out  in  one  step.  The  instrument  red 
noise  part  of  the  measurement  error  is  filtered  by  the  Wiener  smoother  in  the  frequency 
domain.  Thus,  the  gradient  data  needs  to  be  on  a  two-dimensional  grid.  The  application 
of  this  method  for  the  estimation  of  the  gravity  disturbance  vector  at  the  surface  of  the 
earth  from  gradient  data  at  survey  altitude  does  not  provide  accurate  results  at  the  edges 
of  the  area  due  to  spectral  leakage. 

Another  method  that  is  based  completely  on  frequency  domain  approach  was  proposed 
by  Vassiliou  (1985,  1986).  The  method  is  based  on  the  application  of  multiple  input-single 
output  filtering  equations,  using  as  inputs  the  linearly  correlated  second-order  gravity 
gradients  and  as  output  the  first-order  gradients.  In  this  way,  each  first-order  gradient 
is  estimated  from  a  combination  of  its  gradients  in  the  frequency  domain.  The  method 
uses  all  the  gradient  measurements  at  once  for  the  whole  area.  To  make  the  method 
computationally  efficient,  the  Feist  Fourier  Transform  (FFT)  is  employed.  In  this  method, 
gradient  plane  integration  and  downward  continuation  are  performed  in  one  step  using 
FFT  techniques.  This  method,  as  in  any  frequency  domain  technique,  suffers  from  spectral 
leakage.  Additionally,  it  requires  the  data  points  to  be  on  a  regular  two-dimensional  grid 
and  assumes  flat  earth  approximation. 

Parametric  least  squares  estimation  is  the  key  to  another  method  where  the  measured 
gravity  gradients  are  represented  as  a  set  of  parameters  pre-multiplied  by  a  design  matrix 
observed  in  additive  noise  ( Peacock ,  1985;  Center  and  Peacock ,  1985).  The  parameters  can 
be  point  masses,  spherical  harmonics,  two-dimensional  Fourier  coefficients  or  Hankel  trans¬ 
form  coefficients.  In  practice,  the  parameters  used  are  Fourier  series  coefficients.  Thus, 
the  design  matrix  contains  sampled  values  of  the  Fourier  series  basis  functions  evaluated 
at  the  measurement  points.  Singular  value  decomposition  (Lawson  and  Hanson ,  1974)  is 
used  to  decompose  the  design  matrix  into  a  product  of  a  left  orthogonal  matrix,  a  middle 
diagonal  matrix  and  a  right  orthogonal  matrix.  The  elements  of  the  diagonal  matrix  are 
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the  singular  values  of  the  design  matrix.  For  large  sets  of  measurements  a  single  singular 
value  decomposition  is  computationally  impractical.  Therefore,  a  sequential  singular  value 
is  applied.  With  this  sequential  singular  value  decomposition  algorithm,  measurements 
are  processed  in  a  sequence  of  relatively  small  batches.  At  each  stage,  previous  measure¬ 
ment  and  prior  probability  distributions  are  compressed  into  an  equivalent  measurement 
vector.  After  processing  the  new  set  of  measurements,  the  parameter  space  is  orthogo¬ 
nally  decomposed  into  three  sets  of  linear  combinations:  those  that  are  well  estimated, 
those  that  are  partially  estimated  and  those  that  are  essentially  unobserved  based  on  the 
magnitudes  of  the  signal-to-noise  ratios.  When  this  sequential  operation  is  complete,  the 
optimal  estimates  of  the  sequential  singular  value  decomposition  states  are  formed  and  the 
states  can  be  used  to  form  optimal  estimates  of  the  disturbance  vector  at  selected  points. 
The  application  of  sequential  singular  value  decomposition  makes  the  solution  of  least 
squares  problems  with  large  data  sets  computationally  efficient.  However,  to  implement 
a  parametric  least-squares  algorithm,  a  set  of  basis  functions  must  be  explicitly  chosen. 
The  parameters  then  become  the  coefficients  of  a  functional  expansion  in  terms  of  these 
basis  functions.  Since  the  number  of  computations  grows  as  the  cube  of  the  number  of 
parameters,  the  choice  of  parameters  becomes  the  critical  design  issue. 

A  data  processing  technique  that  transforms  a  gradiometer  survey  network  into  an 
electrical  network  was  proposed  by  Eckhardt  (1986).  In  this  method  GGSS  survey  data 
is  analyzed  by  analyzing  an  electrical  network  that  is  isomorphic  to  the  survey  network. 
The  integrated  gradients  between  the  nodes  where  the  survey  lines  cross  correspond  to 
the  applied  voltages  between  the  nodes  of  the  network:  the  gradient  variances  correspond 
to  the  internodal  resistances;  the  elements  of  the  adjusted  gravity  vector  correspond  to 
the  nodal  voltages;  and  the  solution  variances  correspond  to  the  resistances  to  ground. 
Solving  the  electrical  network  is  then  equivalent  to  making  a  least  squares  adjustment 
of  the  survey  network.  The  initial  step  is  to  process  the  raw  GGSS  signals  to  extract 
the  elements  of  the  gravity  gradient  tensor;  next  ground-truth  measurements  are  upward 
continued  to  the  network  tie  points;  then  the  internodal  adjustment  is  done;  finally  the 
gravity  vector  is  calculated  along  the  tracks  between  the  nodes  and  the  field  is  interpolated 
and  downward  continued  to  the  reference  surface.  This  method  differs  from  most  other 
techniques  proposed  in  that  it  is  not  based  on  the  least  squares  collocation  and  thus  does 
not  use  expected  correlations  of  the  gravity  field  and  measurements. 

1.6  Technical  Approach 

Several  different  methods  of  GGSS  data  processing  was  discussed.  Clearly  the  straightfor¬ 
ward  application  of  least  squares  collocation  is  not  a  practical  proposition.  Therefore,  all 
the  methods  approach  the  problem  with  the  idea  of  making  the  computations  tractable. 
The  methods  are  then  not  optimal.  Some  of  the  methods  do  not  use  any  gravity  field 
model  for  the  survey  region  and  those  that  do  are  restricted  to  homogeneous  and  isotropic 
covariances.  But  as  Rummel  and  Schwarz  (1977)  has  shown  the  homogeneous  model  may 


not  always  be  satisfactory  and  the  structure  of  the  anomalous  potential  is  more  realistically 
represented  by  non- homogeneous  weighting  functions.  Morrison  (1975,  1977)  has  shown 
that  the  assumption  of  isotropy  is  not  always  valid  and  Kearsley  (May  1977,  July  1977)  has 
investigated  non-stationary  estimation  where  he  shows  the  superiority  of  two-dimensional 
non-isotropic  covariances. 

Another  fundamental  assumption  in  the  frequency  domain  techniques  is  the  require¬ 
ment  of  stationarity  for  the  gravity  quantities  and  for  the  measurement  noise.  But  as  Nash 
and  Jordan  (1978)  pointed  out,  “Gravity  anomalies  are  neither  stationary  nor  isotropic; 
some  areas  are  rough,  other  smooth,  and  most  areas  contain  linear  (non-isotropic)  fea¬ 
tures  corresponding  to  a  mountain  range,  continental  shelf,  fault,  etc.  The  non-stationary, 
non-isotropic  character  of  gravity  anomalies  is  hardly  surprising  since  approximately  90 
percent  of  the  energy  in  gravity  anomalies  is  caused  by  terrain  (and  the  associated  isostatic 
compensation  at  the  Mohorovicic  discontinuity),  and  most  terrain  contains  obvious  linear 
features.”  Thus  the  frequency  domain  approach,  through  computationally  attractive  really 
does  not  address  the  nonstationarity  and  nonisotropy  of  local  gravity  modeling. 

Nash  and  Jordan  advocate  the  application  of  random  process  theory  to  the  problems 
of  geodesy  particularly  in  the  modeling  and  estimation  of  the  fine  structure  of  the  earth’s 
gravity  field.  They  state  that:  “The  theoretical  basis  for  statistical  models  of  gravity 
anomalies  is  somewhat  tenuous.  ‘Ensembles’  and  ‘probability  distributions’  are  difficult 
to  define,  since  there  is  only  one  earth,  and  gravity  at  a  particular  point  is  a  fixed,  deter¬ 
ministic  quantity  (excluding  transient  effects  such  as  tides)  ....  Nevertheless,  statistical 
geodesy  has  survived  (flourished!)  because  of  the  powerful  conceptual,  mathematical,  and 
computational  tools  afforded  by  the  statistical  approach.”  Random  process  models  of 
gravity  were  first  used  by  Hirvonen  (1956,  1962)  and  Kaula  (1957,  1959)  in  the  late  1950’s. 
These  models  take  the  form  of  autocovariance  functions  (acf’s)  and  spherical  harmonic 
power  spectral  densities  (psd’s).  Statistical  models  for  the  deflections  were  proposed  by 
Levine  and  Gelb  (1969)  without  considering  the  associated  anomaly  models.  Shaw  et. 
al.  (1969)  recognized  that  the  Vening-Meinesz  formulas  provide  a  constraint  between  the 
anomaly  and  deflections;  they  used  this  constraint  to  derive  statistical  models  for  the  de¬ 
flections  from  theoretical  anomaly  models.  In  Particular,  Shaw  and  his  coworkers  proposed 
the  ‘exponential  anomaly  model’  and  the  ‘Bessel  anomaly  model’  and  determined  the  as¬ 
sociated  deflection  models.  In  a  similar  paper,  Kasper  (1971)  proposed  a  ‘second-order 
Markov  anomaly  model’  and  determined  the  associated  deflection  models.  The  models 
proposed  by  Shaw  et  al.  and  Kasper  Eire  ‘self-consistent’,  since  the  anomaly  and  deflection 
statistics  satisfy  compatibility  conditions  that  are  based  on  the  Vening-Meinesz  equations. 
However  these  models  were  somewhat  incomplete  since  even  though  the  constraints  be¬ 
tween  the  anomaly  Emd  deflections  were  recognized,  but  the  constraints  relating  to  the 
undulation  were  ignored.  Jordan  (1972)  took  into  account  these  constraints  and  proposed 
the  ‘third-order  Markov  undulation  model’  and  derived  ‘necessary  conditions’  that  must 
be  satisfied  in  order  for  the  anomaly,  deflection,  and  undulation  correlation  functions  to  be 
physically  reasonable.  Other  self-consistent  models  have  also  been  proposed  (  Tscherning 


and  Rapp,  1974;  Tscherning,  1976).  Some  of  these  self-consistent  models  are  also  mathe¬ 
matically  convenient  insofar  as  the  covariance  functions  aloft  can  be  expressed  analytically 
(Bellaire,  1971,  1972,  1977;  Heller  and  Jordan,  1979). 

Self-consistent  models  are  usually  more  complex  because  they  involve  two  or  more 
gravimetric  uncertainties  (vertical  deflections,  gravity  anomaly,  undulation)  as  opposed  to 
ad  hoc  models  involving  only  one  gravimetric  uncertainty.  The  complexity  of  the  self- 
consistent  models  has  a  blessing  in  that  it  makes  possible  to  infer  a  statistical  model  for 
one  gravimetric  uncertainty  (e.g.  deflections)  from  data  pertaining  to  another  gravimetric 
uncertainty  (e.g.  gravity  anomaly).  However,  even  in  the  case  of  self-consistent  models 
the  autocorrelation  function  (acf)  of  one  of  the  gravimetric  uncertainties  is  required  to  be 
hypothesized  and  the  others  derived  via  the  constraining  equations. 

Some  of  these  ad  hoc  or  self-consistent  models  can  be  expressed  in  state-space  form 
where  the  gravimetric  uncertainties  along  a  line  (trajectory)  in  space  can  be  viewed  as  a 
function  of  time,  rather  than  space.  By  definition,  a  ‘state-space’  model  is  one  that  can 
be  expressed  as  a  set  of  linear  ordinary  differential  equations  excited  by  white  noise.  Such 
models  are  convenient  and  useful  for  covariance  simulation  and  Kalman  filtering.  The 
'third-order  Markov  undulation  model’  of  Jordan  (1972)  is  such  a  model.  However,  this  is 
only  possible  because  of  Jordan’s  assumption  of  an  isotropic  acf  for  undulation. 

As  Nash  and  Jordan  (1978)  admit:  “All  of  the  above  mentioned  models  assume  sta- 
tionarity  and  isotropy,  both  of  which  are  unrealistic  assumptions.  One  critic  has  said 
facetiously,  “Statistical  gravity  models  are  nice  but  they  disagree  with  the  data.”  This 
criticism  is  embarrassingly  valid.  Hopefully,  new  models  will  be  developed  in  the  near 
future  that  account  for  nonstationarity  and  nonisotropy.”  As  a  result  they  state  that  the 
first  major  problem  in  statistical  geodesy  as:  “Nonstationary  nonisotropic  models  need  to 
be  developed  to  replace  the  simple  models  used  heretofore.”  They  further  go  on  to  advo¬ 
cate  the  application  of  modem  estimation  and  control  theory  t.o  geodesy  acknowledging 
though  that:  “The  application  of  modem  estimation  theory  (Kalman  filtering)  to  gravi¬ 
metric  uncertainties  is  complicated  by  the  fact  the  the  quantities  of  interest  are  two  or  three 
dimensional  random  processes,  in  the  sense  that  two  or  three  coordinates  (e.g.,  north-east, 
north-east-down)  are  needed  to  specify  them.”  They  further  state  that:  “The  normal  re¬ 
cursive  approach  ...  is  difficult  to  pursue  directly  for  geodetic  uncertainties  because  they 
are  not  time  series,  but  rather  they  are  two  and  three  dimensional  spatial  random  pro¬ 
cesses.  For  such  processes  it  is  generally  not  possible  to  write  down  linear  differential 
equations  . . .  that  describe  the  dynamic  evolution  of  the  process.  Rather  it  is  necessary 
to  develop  partial  differential  equations,  which  when  excited  by  two-or-three-dimensional 
white  noise,  ‘produce’  the  desired  statistics.  This  is  a  difficult  theoretical  problem;  even  if 
done  correctly,  it  is  not  clear  how  to  solve  the  estimation  problem.” 

Motivated  by  these  remarks  the  approach  taken  here  is  to  exploit  the  marriage  of 
physical  geodesy  and  random  process  theory.  We  shall  solve  Laplace’s  equation  with  the 
unknown  mass  distribution  below  the  surface  of  the  earth  modelled  as  a  two-dimensional 
white  noise  layer  representing  the  vertical  derivative  of  the  disturbance  potential  to  any 
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pre-specified  order.  This  results  in  a  series  solution  of  the  disturbance  potential  wherein 
the  unknown  coefficients  of  the  expansion  are  forced  to  be  uncorrelated  by  invoking  the 
Karhunen-Loeve  condition.  It  is  shown  that  the  disturbance  potential  covariance  obtained 
from  this  model  is  both  non-stationary  and  non-isotropic. 

Next,  the  six  (6)  gravity  gradient  measurements  are  represented  in  terms  of  the  Karhunen- 
Loeve  series  expansion  of  the  disturbance  potential  resulting  in  six  basis  functions.  These 
basis  functions  are  shown  to  be  orthogonal.  A  linear  mean  square  estimator  utilizing  all  the 
gravity  gradient  measurements  simultaneously  is  obtained  in  continuous  domain  by  solv¬ 
ing  an  integral  equation  involving  the  estimator  gains  which  are  represented  by  the  same 
orthogonal  basis  functions.  The  discrete  implementation  of  the  estimator  is  facilitated 
by  exploiting  certain  orthogonality  relationships  of  the  transformation  matrices  such  that 
matrix  inversions  are  not  necessary  at  all.  The  gravity  disturbance  vector  is  obtained  in  a 
two-dimensional  grid  which  can  be  denser  than  the  measurement  grid  and  also  at  any  al¬ 
titude,  including  the  surface  of  the  earth.  Thus,  interpolation  and  downward  continuation 
are  performed  automatically. 

1.7  Overview  of  Report 

This  report  is  organized  as  follows.  Chapter  2  presents  the  gravity  model,  whereas  Chapter 
3  discusses  the  gravity  measurement  model.  The  estimation  algorithm  is  derived  in  contin¬ 
uous  domain  in  Chapter  4.  Discrete  implementation  of  this  continuous  domain  estimator 
is  shown  in  Chapter  5.  Chapter  6  concludes  with  a  brief  summary  of  research  performed, 
highlights  of  research  achievements  and  suggestions  for  future  research. 
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2  GRAVITY  MODEL 

2.1  Introduction  &  Summary 

In  this  chapter  a  local  gravity  model  for  the  survey  region  is  developed.  The  basis  of  the 
model  is  to  represent  the  unknown  mass  distribution  below  the  surface  of  the  earth  as  a 
two-dimensional  white  noise  layer  representing  the  vertical  derivative  of  the  disturbance 
potential  to  any  pre-specified  order.  After  defining  a  rectangular  survey  region,  the  gravity 
field  in  this  region  for  this  local  white  noise  layer  model  is  shown  to  be  a  linear  superposition 
of  a  particular  solution  and  a  complementary  solution.  The  particular  solution  involves  the 
solution  of  Laplace’s  equation  with  zero  boundary  conditions  but  non- zero  sources  specified 
as  a  two-dimensional  white  noise  layer  below  the  surface  of  the  earth.  The  complementary 
solution,  on  the  other  hand,  involves  solution  of  Laplace’s  equation  with  zero  sources  but 
non-zero  boundary  conditions.  The  contribution  from  the  complementary  solution  is  zero 
if  the  boundary  conditions  are  zero.  Only  the  particular  solution  is  addressed  since  it  is 
assumed  that  the  boundary  conditions  are  zero.  The  particular  solution  of  the  disturbance 
potential  is  represented  as  a  series  solution,  where  in  the  horizontal  coordinates  the  solution 
involves  sine  functions  due  to  the  zero  boundary  conditions  and  in  the  vertical  coordinate 
it  has  the  form  of  an  exponential  function  due  to  the  disturbance  potential  vanishing  at 
infinity.  The  unknown  coefficients  of  the  series  expansion  represent  the  random  variation, 
independent  of  any  spatial  variation.  The  condition  of  uncorrelatedness  of  the  random 
variables  is  invoked  thereby  making  the  series  representation  of  the  disturbance  potential 
a  Karhunen-Loeve  series.  An  expression  for  the  autocorrelation  of  the  Karhunen-Loeve 
coefficients  is  obtained  from  the  expression  of  the  white  noise  layer  model  representing  the 
vertical  derivative  of  the  disturbance  potential.  Finally,  an  expression  for  the  covariance  of 
the  disturbance  potential  is  obtained  which  shows  that  the  disturbance  potential  covariance 
obtained  from  this  model  is  both  non-stationary  and  non-isotropic. 

2.2  Local  Region  Spatial  Domain 

A  local  region  of  interest  is  shown  in  Figure  2  where  x  and  y  coordinates  are  on  the  surface 
of  the  earth  and  the  z  coordinate  is  vertical  above  the  surface  of  the  earth.  The  local 
region  of  interest  is  specified  by  the  spatial  domain  D  given  as 


D(x,y)  =  {x,  y;  0<x<A,  0 <y<B) 


(1) 


where 

A  -  length  of  the  survey  region 
B  -  width  of  the  survey  region 

The  survey  measurements  are  taken  at  an  altitude  H  above  the  surface  of  the  earth. 
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2.3 


Local  White  Noise  Layer  Model 

The  gravity  field  in  this  local  region  is  completely  specified  if  the  following  are  known: 

1)  The  boundary  conditions  of  this  region  in  terms  of  the  potential  values  on  the 
boundaries  of  this  region,  and 

2)  The  mass  distribution  density  below  the  surface  of  the  earth  within  this  region. 
The  model  proposed  here  assumes  that  the  boundary  conditions  are  either  known  or 

can  be  estimated,  i.e. 


T(x,y  =  0  ,z)  =  ft(x) 

(2) 

T(x,  y  =  B,z)  =  f2(x) 

(3) 

T(x  =  0  ,y,z)  =  gi(y) 

(4) 

T{x  =  A,y,  z)  =  g2(y) 

(5) 

T(x,y,z  =  ±oo)  =  0 

(6) 

where  T(x,y,z )  is  the  anomalous  potential  or  disturbance  potential 

With  respect  to  the  unknown  mass  distribution  below  the  surface  of  the  earth,  the 
model  proposed  here  is  that  of  a  two-dimensional  white  noise  layer  representing  the  kth 
vertical  derivative  of  the  disturbance  potential  as  shown  in  Figure  2  such  that 
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l&T, 

\dzk^Xuyu 


dkT  } 

Z)  |®z=-D  (*2>  J/2 »  z )  |®z=-D  /  =  —  *2)<Kt/l  —  1/2)  (7) 


0  <  ii,x2  <  i4  0<yuy2<B  z  =  -D 
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where 


(T02  -  variance  of  white  noise  layer 
k  -  order  of  vertical  derivative  of  T 
D  -  depth  of  white  noise  layer 


2.4  Linear  Superposition  Solution 

The  solution  of  the  disturbance  potential  with  boundary  conditions  specified  by  equa¬ 
tions  (2)  through  (6)  and  mass  distribution  modelled  as  in  equation  (7)  can  be  broken 
up  into  two  parts:  the  particular  solution  and  the  complementary  solution  ( Boyc.e  and 
DiPrima ,  1969). 

The  particular  solution  involves  the  solution  of  the  three-dimensional  Laplace’s  equation 
satisfied  by  the  disturbance  potential: 


\dx2  +  dy2  +  dz2 ) 


T  =  0 


(8) 


with  zero  boundary  conditions 


T(x  =  0,y,z)  =  0 


(9) 


T(x  =  A,y,z)  =  0  (10) 

T(x,y=0,z)  =  0  (11) 


T(x,y  =  B,z)  =  0  (12) 

T(x,y,z  =  ±oo)  =  0  (13) 
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Figure  2:  Local  White  Noise  Layer  Model 
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but  non- zero  sources  specified  as  a  two-dimensional  white  noise  layer  such  as 


'  dkT 

dzk 


(XllVll  2)  |  <j)Z: 


dkT 

=  -D  X-^j(x2,y2,z)  |@z= 


d  f  =  <rl6(xi  -  x2)£(yi  -  y2)  (14) 


0  <  ii,:r2  <  .A  0  <  j/1,1/2  <  B  z  =  —D 


The  complementary  solution  involves  the  solution  of  the  three-dimensional  Laplace’s  equa¬ 
tion  satisfied  by  the  disturbance  potential: 


[ dx 2  +  dy>  +  dz 2  J 


T  =  0 


(15) 


with  zero  sources 


T  =  0 


(16) 


0<j<i4  0  <  y  <  B  —00  <  z  <  +00 


but  non-zero  boundary  conditions 


T{x,y  =  0  ,z)  =  f\ (x) 

(17) 

T(x,y  =  B,z)  =  /2(z) 

(18) 

T(x  =  0,y,z)  =  yi(y) 

(19) 
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T(x  =  A,y,z)  =  g2(y) 


(20) 


T(x,  y,z  =  ±oo)  =  0 


(21) 


In  the  sequel,  only  the  particular  solution  is  addressed.  The  complementary  solution  is 
zero  if  the  boundary  conditions  are  assumed  zero.  The  rest  of  this  development  assumes 
zero  boundary  conditions  thereby  eliminating  the  need  to  address  the  complementary 
solution  any  further. 

2.5  Series  Representation  of  Potential  Function 

The  basic  approach  to  the  solution  of  the  problem  is  obtained  by  representing  the  scalar 
field  T(x,y,z )  in  a  series  expansion  as  follows 


oo  oo 

T(x,y,z)  =  V  (  )  *  Vi 

m=l  n=l 


(22) 


The  above  representation  of  T(x,y,  z)  implies  that 

M  N 

T(x,  y,  z)  —  l.i.m.  amn<fimn(x,  y,  z) 

M  — »  OO  m  =  l  n=i 

N  oo 


(23) 


where  l.i.m.  is  the  limit  in  the  mean  square  sense  such  that 


l.i.m. 
M  — y  oo 
N  —*  oo 


C  M  N 

E  %  |  T(x,  y,z)  —  y  ~  !/i  z) 

l  m=l  n=l 


(24) 


where  E  is  the  statistical  expectation  operator  and  the  functions  4>mn(x,  y,  z )  are  complete 
and  square  integrable. 

Reflecting  on  equation  (22),  we  observe  that  by  this  representation  we  have  separated 
the  random  variation  from  the  spatial  variation  of  T{x,  y,  z).  The  randomness  of  T(x,  y,  z ) 
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is  included  in  amn  which  is  a  random  variable  but  independent  of  any  spatial  variation. 
The  spatial  variation  of  T(x,  y,  z),  on  the  other  hand,  is  completely  included  in  <f>mn(x,  y,  2) 
which  is  spatially  varying  but  statistically  a  deterministic  function. 

2.6  Particular  Solution  Spatial  Functions 

The  solution  of  equation  (8)  with  zero  boundary  conditions  given  by  equations  (9)  through  (13) 
is  given  by  using  the  representation  of  equation  (22)  as 


T(x,y,z) 


2 

Sab 


EE 


amn  sin  amx  sin  bnye  Cm"|z+D| 


(25) 


where 


mn 


Cmn  =  «  +  *«)’ 


(26) 


(27) 


(28) 


and  amn  are  the  unknown  random  coefficients  of  the  series  solution. 

2.7  Karhunen-Loeve  Uncorrelatedness  of  Series  Coefficients 

The  series  representation  of  the  disturbance  potential  given  by  equation  (25)  is  represented 
as  a  Karhunen-Loeve  Series  ( Davenport  and  Root,  1958)  by  enforcing  the  condition  of 
uncorrelatedness  of  the  random  variables  amn  implying  that  (Papoulis,  1965). 


E(otmnOtkl)  —  ^mn&mk&nl 


(29) 
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where  A  mn  is  the  correlation  of  the  random  variables  am7l  and  the  Kronecker  delta  function 
is 


6ij  ~  {  C  ijt  j  (30) 

2.8  Autocorrelation  of  Random  Variables 

An  expression  for  the  correlation  of  the  random  variables  is  now  obtained  by  utilizing  the 
condition  of  uncorrelatedness  of  the  random  variables  and  the  two  dimensional  white  noise 
layer  source  at  depth  D. 

From  equation  (25) 


dkT  =  2(— sgn  (z  +  D))k 

dzk  Jab 


oo  oo 

5Z  52  amnCkmn  sin( am  j  )  sin ( 6n y ) e ~ c'm"  + 01 

m  —  \  »— 1 


(31) 


Substituting  equation  (31)  in  the  left  hand  side  of  equation  (7)  yields 


( dkT  ) 

E  |  ^  *  Qzk  1/2)  2)  laz=-D  r 


OO  OO  OO 


=  js  52  525252  £(a- 

m=l  n=l  fc=l  /= 1 


»a*/)cmnc*/sin  am*i  sin  bnyi  sinafcx2  sint;y2 


(32) 


Substituting  equation  (29)  in  equation  (32)  reduces  the  quadruple  summation  to  double 
summation  such  that 


a*T, 


dkT, 


dzk-(x \,y\,z)  |o *=-£>  x^-(x2,y2,z)  |o*=-d| 
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(33) 


A  oo  oo 

m  =  l  n=l 


c2k 

mn^rrjn 


sin  amX\  sin  bny i  sin  amx2  sin  bny2 


4 

AB 


Amn4- 


Ik 

mn 


oo  oo 

^2  sin amxi  sinamx2  ]T  si n&„y2 

m=l  n=l 


(34) 


where  A mnc™n  is  assumed  constant.  Now,  using  the  series  representation  of  Dirac  delta 
functions 


9  °° 

-  x2)  =  Yi  sinam3-i  sinamx2  (35) 

A  m  =  1 


and 


2  00 

b(y i  -  Vi)  -  5Zsin6nyi  sin6ny2 


(36) 


in  equation  (34)  yields 


( dkT  dkT  1 

E  {  0^(*l,yi,z)  loi=-D  *-gZk(x  2,V2,z)  |®2=-D  ? 


=  \mnC%Jixt  -  x2)S(y ,  -  y2) 


(37) 


Substituting  equation  (37)  in  equation  (7)  gives 


KnC^Jix i  -  x2)S(yl  -  y2)  =  cr£6(xt  ~  *2)%!  -  y2 ) 


(38) 


which  leads  to 


A 


mn  — 


(39) 
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This  is  a  consequence  of  assuming  a  white  noise  layer.  In  principle,  this  assumption  could 
be  relaxed  to  permit  any  covariance  representable  in  the  form  given  by  equation  (33).  It 
should  be  obvious,  however,  that  bounded  stationary  covariances  are  ruled  out  by  the  zero 
boundary  condition. 

2.9  Non-isotropic,  Non-stationary  Covariance 

An  expression  for  the  covariance  of  the  disturbance  potential  is  easily  obtained  by  using 
equations  (25)  and  (29)  as  follows: 

Rtt{xi  ,  yx ,  zi ;  x2 ,  y2 , 22 ) 


=  E  {T{xuyuZi)T(x2,y2,z2)} 


f  2  00  °° 

=  E  \  ~7=f=  Y1  X  amn  sin  “ml  1  sin  6ny1e'c"'"l*1+D| 

l  v  AB  m=1  n=l 


r)  OO  OO 

x  X  X  Qu  sin  a*;r2 sin  biy2e'Ck‘lZi+D{ 
vAB  fc=1 


^  OO  OOOOOO 

=  Tp  X  X  XX£(a"*na*/)sinamr,  sin&nt/i  sin  akx2  sin  fyy2e~c,n"|ei+D|~ct‘'Z2+D| 

Ati  m  =  l  ij=1  1=1 


4  OO  OO 

=  51  XI  Amnsinamx,  sin  6nyi  sinamx2  sin  &ny2e~Cm"(l*1+D|+l*2+r,1)  (40) 

m=l  n= 1 


0  <  xi , xj  <  A  0  <  yi,y2  <  B  -  D  <  zuz2  <  +00 

Note  that  for  zi,z2  >  —D  the  vertical  variation  is  given  as  (z,  +  z2  +  2D)  which  has  the 
right  form  as  shown  by  Moritz  (1980)  [page  171]. 


3  SENSOR  MODEL 


3.1  Introduction  &  Summary 

In  this  chapter  the  mathematical  model  of  the  gradiometer  sensor  is  examined.  Beginning 
with  the  basic  definition  of  the  three  (3)  inline  and  three  (3)  crossline  gradiometer  out¬ 
puts  in  terms  of  the  disturbance  potential  a  set  of  six  (6)  basis  functions  are  defined  in 
terms  of  the  Karhunen-Loeve  series  expansion  of  the  disturbance  potential.  These  basis 
functions  are  shown  to  be  orthogonal  within  the  survey  region  both  individually  as  well 
as  cumulatively.  An  expression  for  the  autocorrelation  function  of  the  six  (6)  gradiometer 
signals  taken  together  is  obtained  in  terms  of  the  basis  functions  and  the  autocorrelation 
of  the  Karhunen-Loeve  coefficients.  This  is  then  used  to  derive  an  integral  equation  for  the 
gradiometer  signal  autocorrelation  function.  The  measurements  of  the  inline  and  crossline 
gradients  are  contaminated  with  additive  noise  whose  autocorrelation  function  is  assumed 
to  be  a  diagonal  matrix. 

3.2  Gradiometer  Signal  Representation 

The  gradiometer  sensor  signal  S(x,y,z)  has  six  (6)  components:  three  (3)  inline  signals 
and  three  crossline  signals.  The  six  (6)  signals  of  the  gradiometer  are  given  below 


Si(x,y,z) 


S2(x,y,z) 


5j(x.  9,  *) 


1  fd2T  _  d'2T\ 

2  \  dx 2  dy 2  / 


1  (d2T  d2T\ 

2  [dy*  ~  dz* ) 


1  (dn  _  d2T\ 

2  \  dz 2  dx 2  ) 


(41) 


(42) 


(43) 


SA{x,y,z) 


d2T 

dxdy 


(44) 


Ss{x,y,z) 


d2T 

dydz 


(45) 
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Se(x,y,z) 


d2T 

dzdx 


(46) 


3.3  Definition  of  basis  functions 

Taking  the  appropriate  partial  derivatives  of  equation  (25)  and  substituting  in  equa¬ 
tions  (41)  -  (46)  gives 


OO  OO  (  _ J 

Si(x,y,z)  =  53  X3  <*mn  sin  amx  sin  bny  |^^(am  ~  bl) 


}"Cmn  \Z^~D\ 


OO  OO  <  1 

^2{x,y,z)  =  53  5  1  Qmn  sin  amx  sin  bny  (cmn  -I-  frn)e  ^  j 


oo  oo  r  j 

S2(x,y,z)  =  53  Xj<*m«»inom*sin6kly|-^=(c2IfI  +  <£)e 


2  U'Cmn  |2+ 


*?4  y,  2)  —  ^  ^  Q:mn  COS  OitjX  COS 


b  v  < —2=a  b  e~Cmn,z+Di  \ 

nV\jABmn  j 


S5(x,y,z )  =  53  53amnsinamxcos6ny|-^^L£l6nCmne  Cm"l*+D'|  (51) 


OO  OO  ( 

S&(x,y,z )  =  53  J2amnCOsa”'xsmbny  i 

m= ]  n=l  l 


—2  sgn(z  +  D)  Cm„|z+D| 

VlB  mn 
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a t  z  =  H  these  become 


I  -i 


4 


S\(x,y)  =  ^  ^amnsinama;sin6ny  |-^==(a^,  -  6^)e  Cmn|//+D[| 


(53) 


OO  OO  f  _  1  'I 

S2(x,y)  =  53  51  atmnsmamxsrnbny  jyj^(cmn  +  6^)e-c"'"|H+D||  (54) 


oo  oo  C  ^ 

■Ss^y)  =  Y1  ]C<imnsinama:sin6ny  |-^=(c^n  +  a^)e 


2  '|c-Cmn|W+D|l 


'VmnJ 


OO  OO  (2  1 

S4(x,y)  =  53  Xjamncosamxc°s6ny  j  y^am&"e_Cmn|H+  J 


'Ymn4 


(55) 


(56) 


s5(z,  y)  =  53  5>musinamzcos 

m=  1  n=l 


Ky  {  ~2-8|^+  P)-inc„e—  (57) 


OO  OO  f 

s6(x,y)  =  53  5Za™"cosa"'xsin^y  i 

m= I  n= 1  V 


— 2sgn  {H  +  D)  - cm„\H+D | 

umLmnc 


v/lB 


(58) 


Representing  ail  six  (6)  signals  of  equations  (53)  through  (58)  as  a  vector  signal  S(x,  y) 
such  that 


S(x,y)  —  [Si(*,y)  S2(x,y)  S3{x,y)  S4{x,y)  S5(x,y)  56(z,y)]r 


(59) 
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we  obtain 


OO  OO  6 

S{x,y)  =  52  S  52  CXmn1pmni{xi  i/) 

m=l  n=l  *=1 


(CO) 


OO  OO 

—  )  ^  y  ^  Qmn^mw(^i  2/) 

m=l  n=l 


‘ 


where  the  individual  basis  functions  t/wi  are  defined  as 


(61) 


lT 


V^mnl 


sina,„xsm 


Ky  | 


-l 

Tab 


(<£  -  6^)e-Cm"i"+°i  i  ooooo 


(62) 


V’mnZ 


0  sin  amx  sin  bny  | 


+  &2)e-c""0H+D|  l  0  0  0  0 


(63) 


t^mn3 


0  0  sinamzsin&„y  +  am)e  C"’"|H+D|j  000 


(64) 


0mn4  — 


0  0  0  cosamxcosfe, 


'ny  | 


— ■'*"}  00 


(65) 
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1pmn5  — 


00  00  sin  am  x  cos  bn  y  j  — 2  S  ^  bncmne  C”*»IH+Dlj  0 


(66) 


^mn6  — 


0  0  0  0  0  cos  amx  sin  bny  j 


-2sgn(ff +  £>) 
y/AB 


-cmn|tf+D|  l 


(67) 


and  the  cumulative  basis  function  ipmn  is  the  sum  of  the  individual  basis  functions  4'mni 

e 

4'mn{x,  y)  =  Yl  4>mni(r,  y)  (68) 

1=1 

3.4  Orthogonality  of  Basis  Functions 

The  particular  choice  of  basis  functions  given  by  equations  (62)  -  (67)  is  particularly 
desirable  because  they  are  orthogonal.  Orthogonality  of  these  basis  functions  is  examined 
by  evaluating  the  integral 


/  4>kij(x>y)4’mni(x,y)dxdy 
Jy=0 


for  each  combination  of  i,j  where  i  ~  1,2, 3, 4, 5, 6  and  j  =  1,2, 3, 4, 5, 6.  For  example, 
consider  j  =  1  and  i  =  2,  then 


Jrx=A  ry=B 

/  [t*/i  sin  OkX  sin  b/y  0  0  0  0  0] 

x=0  J  y=0 


o 

Imni  sin  amx  sin  b„y 
0 
0 
0 
0 


dxdy 


=  0 


(69) 
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But,  when  t  =  j  =  1,  then 


L 


v=B  b 

sin  bny  sin  biydy  =  —6nt 
y=o  L 


(72) 


cos  amx  cos  atxdx  =  — 


(73) 


L 


y=B  B 

cos  bny  cos  btydy  =  -~6n, 

y=0  i 


(74) 


where  am,ak,bn,bi  are  defined  in  equations  (26)  -  (27).  Thus,  individual  basis  functions 
defined  in  equations  (62)  through  (67)  are  easily  seen  to  be  orthogonal  satisfying  the 
orthogonality  condition  given  by 


rx~A  ry—B 
J  X—0  •/ yrrQ 


AB 

li(x,y)dxdy  =  'Ykljlmnibmkbnlbij 


(75) 


The  orthogonality  of  the  cumulative  basis  functions  is  examined  as  shown  below. 


[Z~A  ry=B 

I  n  /  Vmn^kidxdy 

r=0  >7  y=0 


rx=A  ry=B  .. 

—  /  [7mm  smamxsin&„y  7mn2  sin  amx  sin  bny  7mn3  sin  amx  sin  bny 

J  X=Q  J  y=0 


^fmn 4  COS  OmX  COS  bny 


”Tmn5  sill  Grn  %  COS  6ny 


Ttnn6  cos  sin  6ny] 


X 


7fcii  sin  a* i  sin  b\y 
Ikn  sina/txsin  6jy 
7k<3  sin  afcx  sin  6;y 
~fkiA  cos  a*  i  cos  b;y 
jkis  sin  akx  cos  i/y 
7*n  cos  a*  i  sin  biy 


dxdy 
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=  7mni7 ux  r=A  sin  amx  sin  akxdx  f  sin  bny  sin  b,ydy 

Jz—O  _  ^  *  u  _ _  ✓ 

fim*  fi»* 


+7mn27fci2  /  sin am X sin aixdr  /  smbnysmbiydy 

J  x=0  _  ^yscil  _ 


f  irn* 


+lmnslkiz  J  q  sin amx  sin  akxdx  smbny  smb,ydy 

fi-n*  f*»< 


+7mn47fci4  J  q  cos  ctmx  cos  o.kxdx  cos  6„y  cos  f>fydy 


7S" 


fin. 


+7mn57fc/5  J  o  sin  amxsin  akxdx  cos  b„y  cos  b,ydy 


f«n. 


+7mn67fci6  /  cos  amx  cos  Ofcxdi  J  sin  6nysm  bjydy 
7i=o  _ _ ^ 


fin. 


a 

=  7mnl7fcn-J-0mfc0„|  +  7mn27*<2  4  «mfcOn/ 


.  ,  AB 

+7mn37fci3-J-^m*0„|  +  ImnilklA  ^ 
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AB  f  AB .  x 

“t “]mr,5^kl5  “  OmkOnl  't'  7mn67M6  .  Omk*n 

4  4 


^  {^mnl  d”  7mn2  d~  7mn3  d~  Tmn4  d”  %nn5  d”  ^mn6 1  ^mk^nl 


.  ^  'Ymni~tkU&mk&nl 

4  i=l 


where  the  7’s  jure  defined  in  equations  (62)  -  (67).  Thus,  the  cumulative  basis  functions 
are  also  orthogonal  satisfying  the  orthogonality  condition  given  by  equation  (76). 

3.5  Measured  Signal  Autocorrelation  Function 

An  expression  for  the  autocorrelation  function  of  the  gradiometer  signals  can  be  obtained 
using  equations  (60)  and  (29)  as  shown  below: 


•ffss(*l,J/i;*2,y2) 


=  E  |5(xj,  t/i  )ST(i2,  j/2)} 


00  00  6 


=  £  E  ]££amnlk»n.-(*t,Vl)  ]£  £  £  Qkl^L}(X2,  Vi) 


a=I  n=l  i  =  l 


*=1 /=■! j= 1 


00  00  6  00  00  6 


=  E  E  E  E  E 

m= 1  n=l  i=l  ik=l  /=1  j= 1  ~  ^ 


(77) 


oo  oo  6  6 

=  2D  2D  2D  2D  *mn*mni(*l,Vl)*ZnA*2,V2) 

m=l  n=l  i=l  j—l 

where  Amn  was  defined  in  equation  (39)  and  t/>mnl  is  defined  in  equations  (62)  -  (67). 

3.6  Signal  Autocorrelation  Integral  Equation 

An  integral  equation  for  the  measured  signal  autocorrelation  function  R$s  can  be  obtained 
using  equations  (77)  and  (75)  as  shown  below: 


/x2=/t  r\n-B 

/  /  Rss{Xl,yi',X2,y2)i>mni(X2,y2)dx2dy2 

J  X2  =0  J  =0 


/X-2=A  ry~2~S  ^  ^  ^  ^ 

/  2D  2D  2D  2D  ,  yi)V>*/>(*2,  y2)dx2dy2 

-2=0  JW=°  fc=U=l£=lj=l 


6  ®  rzi=A  rV2~t> 

~  E2D2DE  Xkrtkit(xi,yx)  /  /  ipki}(x2,y2)il>mn,(x2,  y2)dx2dy2 

k=i  1=1  <=1  ;sl  •/l2=°  •/»=° 


rj=4  rtn~B 


oo  no  6  6 

/=!  t=l  j-1  4 


=  -7-Amn7^n,5DV’mn£(a:i,yi)  (78) 

4  i=i 

3.7  Measurement  Noise  Representation 

The  measurement  Z(x,y)  at  z  =  H  of  the  gradiometer  signal  S(x,y)  is  contaminated  with 
an  additive  noise  V(x,y)  such  that 
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Z(x,y)  =  S(x,  y)  +  V(x,y) 


(79) 


The  autocorrelation  function  of  the  additive  noise  is  assumed  to  be  a  diagonal  matrix  such 
that 

Rvv  (xuyi\x2,y2)  =  {V(xx,yl)VT(x2,y2)}  (80) 


0  0  0  0  O' 

0  o\  0  0  0  0 

0  0  <t?  0  0  0  c.  sc,  .  .... 

0  0  0  ^00  <5(3rj  -  x2)6(yi  -  y2)  (81) 

0  0  0  0  a\  0 

0  0  0  0  0 

The  vector  signal  S(x,  y)  and  measurement  noise  V(x,  y)  are  assumed  to  be  uncorrelated, 
i.e. 

£{5(x,,yi)FT(x2ly2)}  =  {V(x1,y,)S/’(x2,y2)|  =  0  (82) 
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4  ESTIMATION  ALGORITHM 

4.1  Introduction  &  Summary 

In  this  chapter  an  estimation  algorithm  is  developed  in  the  continuous  domain.  Thus, 
herein  it  is  assumed  that  measurements  are  available  in  continuous  form  throughout  the 
survey  region.  A  linear  mean  square  estimation  is  proposed  from  which  it  is  shown  that 
estimating  the  disturbance  potential  or  any  functional  of  the  disturbance  potential  such  as 
the  gravity  vector  or  the  gravity  gradients  is  equivalent  to  simply  estimating  the  Karhunen- 
Leove  coefficients.  Exploiting  the  orthogonality  principle  of  linear  mean  square  estimation 
and  the  uncorrelatedness  of  measurement  noise  and  signal  leads  to  a  vector  integral  equa¬ 
tion  involving  the  matrix  autocorrelation  function  of  the  signals  to  be  estimated.  This 
vector  integral  equation  involves  the  estimator  gains  which  are  represented  by  means  of 
the  same  orthogonal  basis  functions  with  the  unknowns  being  the  coefficients  of  the  ex¬ 
pansions.  An  explicit  solution  of  the  vector  integral  equation  is  then  obtained  which  leads 
to  a  closed  form  solution  of  the  coefficients  in  the  estimator  gains.  Thus,  the  estimates  of 
the  Karhunen-Loeve  coefficients  are  expressed  in  terms  of  the  coefficients  of  the  estimator 
gains  and  measurements  integrals  of  the  six  (6)  gravity  gradient  measurements.  A  key  re¬ 
sult  is  that  the  Karhunen-Loeve  coefficients  axe  estimated  utilizing  all  the  gravity  gradient 
measurements  simultaneously.  Once  the  Karhunen-Loeve  coefficients  are  estimated,  the 
disturbance  potential,  gravity  vector  components  or  gravity  gradients  can  be  estimated  by 
utilizing  their  functional  relationships  to  the  Karhunen-Loeve  coefficients. 

4.2  Linear  Mean  Square  Estimation 

The  desired  smoothed  estimate  5(x,y)  for  the  six-dimensional  vector  signal  S(x,  y)  given 
by  equation  (59)  is  obtained  by  minimizing  the  mean  square  error  between  the  signal 
5(x,y)  and  its  estimate  S(x,y)  under  a  performance  index  J(x,y)  defined  as 


J(x,y)  =  E  ||5(x,  y)  -  S(x,y)|2} 


(S3) 


The  optimal  estimate  S(x,  y)  which  minimizes  the  performance  index  defined  by  equa¬ 
tion  (83)  for  every  point  (x,  y)  in  the  domain  D  defined  in  equation  (1)  at  the  measurement 
altitude  z  —  H  is  given  by  ( Papoulis,  1965) 


5(x,y)  =  E  {S(x,y)\Z(u,v)} 


(84) 
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0<uOi  0  <v<B 


where  the  measurements  Z(u,v,H)  are  taken  over  the  entire  domain  D(u,v )  a*  altitude 
H.  From  equation  (60)  we  deduce  that 


OO  OO  6 

m=l  n=l  i=l 


(85) 


from  which  we  conclude  that  obtaining  S(x,y)  is  equivalent  to  simply  estimating  amn  from 
the  given  measurements  Z(u,v).  Tf  we  assume  a  linear  estimator  this  can  be  written  as 


=  rA  rB 

J. r=0  J  y=0 


ZT(x,y)Kmn(x,y)dxdy 


(86) 


where  dmn  is  a  scalar  and  Z(x,  y)  and  Kmn(x,  y)  are  each  6x1  column  vectors  representing 
the  survey  measurements  and  the  estimator  gains  respectively.  In  the  case  of  Gaussian 
statistics,  the  linear  estimator  given  by  equation  (86)  is  equivalent  to  the  general  estimator 
given  by  equation  (84). 

4.3  Application  of  Orthogonality  Principle 

For  linear  mean  square  estimation  the  orthogonality  principle  is  given  by  ( Papoulis ,  1965) 

E^Z{xuyi)  [5(12,2/2)  -  5(x2,2/2)]T|  =  0  (87) 


Substituting  equations  (60)  and  (85)  in  the  left  hand  side  of  equation  (87)  leads  to 


E  |z(ii,yt)  [5(12,2/2)  ~  5(:r2, 2/2)]  } 


E  |z(xi,j/i) 


00  00  6  00  00  6  I 

yi  5"!  y  OtmnlftmniiZi,  2/2  )  ~  ^  '  dmnt/’mnj(X2,  t/2) 

.m=l  n=l  «  =  1  rn=l  n=l  t  =  1  J 
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=  E  <  Z(xi,yx) 


-i  T' 


E  E  E(a™«  - 

m=l  n=l  t=l 


i)lpmni(X2,y2) 


{  oo  oo  6 

^..yOEEEK 

m=l  n=l  i=l 


=  E 


Z(xuyx){ar, 


-ow)^n.(x2>y2) 


-  amn)ip*nt{x2,y2) 


} 

} 


oo  oo  oo 

=  E  E2  El  ■^'{■^(Xl'yi)(fkmn  —  ttmn  )}  xpmni(x2,  2/2  ) 

m=l  n=l  i=l 


—  E  [Z(X\,y\ )(Ofmn  ®mn)} 


which  implies 


E  {Z(xi,yi)am„}  >  E  {^(xi>yi)dmn) 


(88) 


4.4  Uncorrelatedness  of  Measurement  Noise  and  Coefficients 

Substituting  the  signal  representation  given  by  equation  (60)  into  the  left  hand  side  of 
equation  (82)  leads  to 


E{V(Xl,yx)ST(x2,y2)} 


=  E 


V{xi,yx) 


oo  oo  6 

EEE 

,m  =  l  n=  1  1=1 


&mn‘,ftmni{3'2i  1/2  ) 


{OO  OO  6  ^ 

J2Y,amn^lnt{x2,y2)\ 

m=l  n=l  t  =  l  J 


oo  oo  6 

=  Y,  YYE  {V(X^yi)°‘mrl}i'ln,(X2,y2) 

m=l  n=l  i=l 


=  E  {V(xuyt  )amn} 


which  implies 


^{VVi^iKnn}  =  0  (89) 

4.5  Integral  Equation  of  Estimator  Gains 

A  vector  integral  equation  involving  the  matrix  autocorrelation  functions  of  the  vector 
signal  S(x,y)  and  the  vector  estimator  gains  Kmn(x,y)  can  be  obtained  as  follows: 

Using  equations  (79),  (89),  (60)  ,  and  (29),  the  left  hand  side  of  equation  (88)  simplifies 
to 


E  {Z(xx,yx  )a 

mn  } 


=  E  {{S(xi,yt)+  V(xuyi)]amnj 


=  E  {5(ii,yi)amn}  +  E{V(xu  !/,)amn) 


=  E 


oo  oo  6 
k=ll=\j  =  \ 


Q'mn 
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oo  oo  6 

mn 

k=\  1= l  j=i 


oo  00  6 

=  12  Y1  H  ^klf>mkfinlll’klj(Xl ,  S/l  ) 
*=1  /=1  >=1 


6 

=  ^mn  ,  J/l  ) 

J  =  1 


(90) 


Using  equations  (86),  (79),  (82)  ,  and  (80),  the  right  hand  side  of  equation  (88)  simplifies 
to 


E{Z(xx,yi)a 

mn  } 


{rx2=A  ryi  =  B  } 

z(xi,yi)  /  Z {x2,y2)Kmn(x2,  y2)dx2dy2  > 

Jl2—  0  J  )J7  J 


(  [* 2=A  rsn=B 

=  E\  / 

(yi2=0  •'1/2=0 


Z(xuyi  )Zr(x2,  y2)A'mn(x2,  y2)dx2dy2 


} 


E 


f  ri2=A  ry2~B 

/ 

^«/l2=0  ^  1/2=0 


[5(xi,y,)  +  Vr(x,,y,)][5(x2,y2)  +  U(x2,y2)]r  Kmn(  X2,yi)dx2dy2 


~  L  /  =0  yi)-S’XC^2,  y2>}  +  E  {5(xi,y,)UT(x2,y2)}  +  E  {V(x,,  y,  )S‘  (x2,  y2} 
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+  E{F(Il,y1)Vr7’(x2,y2)}]  Kmn(x2,yi)dx2dyi 


{/?5s(xi,yi;x2,y2)  +  ^vv(xi,yi;x2,y2)}  Kmn(x2,y2)dx2dy2 


(91) 


which  is  known  as  the  Wiener-Hopf  equation.  Equating  equations  (90)  and  (91)  yields 


6 

xi,yi;x2,y2)  +  ^vv(xi,yi;x2,y2)}  Kmn(x2,y2)dx2dy2  =  Amn  ][^mnj(*i,yi) 

2  =  1 


(92) 


4.6  Orthogonal  Representation  of  Estimator  Gains 

The  vector  estimator  gains  I\mn(x ,  y)  of  the  linear  estimator  given  by  equation  (86)  can  be 
represented  by  the  orthogonal  basis  functions  given  by  equations  (62)  through  (67)  such 
that  (Morse  and  Feshbach,  1953) 


Am„(*,y)  =EE£/^"iM*,y)  (93) 

jt=i  i=i  j=i 

4.7  Solution  of  Estimator  Gains  Integral  Equation 

The  estimator  gains  Kmn(x,y)  given  by  equation  (93)  have  the  basis  functions  tpktj(x,y) 
which  are  known  as  given  in  equations  (62)  through  (67)  but  the  /?£)"’ s  are  yet  to  be 
determined.  An  explicit  expression  for  the  /?^"’s  is  obtained  by  solving  the  vector  integral 
equation  (92)  as  follows. 

Using  equations  (93)  and  (78)  in  the  first  half  of  the  left  hand  side  of  equation  (92) 
yields 


rii=A  fV2  =  S  , 

/  /  Rss{Xi,y\\X2,yi)Kmn{X2,y2)dx2dy2 

J I2  =0  *'V2=  0 


(T1=A  ryi=B  [~  “  «  1 

=  /  /  Rss(x\,yuX2,y2){2Z22z2^ki}^‘3(x2,y2)}dx2dy2 

Jr,=0  =  0  l  *=,;=!;=!  J 
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'  a\  0  0  0  0  0  O' 

0  a\  0  0  0  0  0 

0  0  03  0  0  0  0  00006 

000**000  EELW.^0 

0  0  0  0  a\  0  0  k=u=ij=i 

0  0  0  0  0  a\  0 

0  0  0  0  0  0  a\ 
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(95) 


oo  oo  6 
jt=i  1=1  j=i 


The  necessity  of  a  diagonal  representation  for  the  measurement  noise  covariance  Rvv 
is  clear  in  the  derivation  of  equation  (95).  Unless,  Ryy  is  diagonal,  the  last  step  of 
the  derivation  of  including  the  variances  a *  within  the  summation  is  not  possible.  This 
ability  to  include  the  variances  within  the  summation  sign  greatly  enhances  the  analytical 
tractability  of  the  estimator  gains.  Substituting  equations  (94)  and  (95)  in  equation  (92) 
leads  to 


AB 

4 


EEEE/^^*«(*i.vi) 


k=l  /=!  ;=1  t=l 


oo  oo  6 


6 


+  E  EE  VjPkiji’kijix  I,  yi )  —  ^mn  E^mn  ;(X!,yi) 

k=ll=lj=l  ^=1 


(96) 


which  can  be  written  as 


CO  OO  [  6 

EE  E 

k= 1  1=1  (j=l 


AB 


Pkij^nltj  E  iM*  1 » yi )  +  ct'WJMxu  yi) 


(=1 


=  An 


lE^- 

;=» 


mnj 


(zi,Vi) 


(97) 

Since  the  right  hand  side  of  equation  (97)  contains  no  summations  on  k  or  /,  a  uniques, 
solution  to  equation  (97)  exists  if  and  only  if  m  =  k  and  n  =  l  such  that 


0k,;  =  /  Wrnfc* 


nl 


_  f  0mnj 

1  o 


when  m  —  k  and  n  =  l 
when  m  jh  k  or  n  l 


(98) 


Using  equation  (98)  in  equation  (93),  the  estimator  gains  A'mn(z,y)  reduce  to 


oo  oo  6 

Am„(*,y)  =  EEE 


k=  1  (=1  j= 1 
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oo  oo  6 

=  y  l  y  y  ftm nj^mk^nl^Pkl}(x  i  J/) 

k= 1  1= 1  j=l 


6 

=  yi  1/) 

i=i 


(99) 


An  expression  for  the  /?mnj’s  in  equation  (99)  is  obtained  by  substituting  equation  (98)  in 
equation  (97)  and  solving  the  resulting  equation  given  by 

6  ( 6  1  6 
^  1  »  7"  ftmnj  •^mn'Ymnj  ^  V,mn((1 1 1  !/l)  “b  ^mnjV’mnj(^l  >  3/l  )  f  =  ^mn  i  Vl  ) 

j  =  l  l  4  (=1  J  i=l 

(100) 

Multiplying  the  second  term  of  the  LHS  by  )C?=i  >  renaming  the  dummy  index  j  to  t  on 

the  RHS,  and  rearranging  yields 

5Z  j]C  {^VtL,  +  VjJ^mnj  -  ^mnj  */Wt(*l ,  Vl )  =  0  (101) 


By  the  linear  independence  of  the  il>mnt(xuy\)  equation  (101)  reduces  to 

6  'AB 


£,  {-^-Amn7mnj  +  Pmn,  =  K 


(102) 


which  can  be  written  in  matrix  form  as 


AB  \  -,2  .  n2  AB\  -v2 

4  ^mnfmnl  '  4  /'mn7mn 2 

AB  \  2  M\  -v2  4- rr2 

4  /Amn7mnl  4  AmnJmn2  •  a2 


L  ^*mn7i 


mnl 


4S  \  -y2 

4  ^mnVrnnS 

AB  \  -,2 

4  Amn7mn6 


4B 


^mnl 

^mn 

Pmnl 

= 

^mn 

L  ^mn6 

^mn 

(103) 
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By  using  the  method  of  Gaussian  elimination,  equation  (103)  can  be  solved  to  yield 


$mnt  — 


(104) 


which  can  be  also  written  as 


0mni  — 


ii?=1^  +  ma 


"'Tnn  (^j  =  l  Tmnjn®=1  %k^j<y£) 


(105) 


where  the  7mnj’s  are  given  in  equations  (62)  -  (67),  Amn  is  defined  in  equation  (39)  and 
the  ays  are  introduced  in  equation  (81). 

4.8  Continuous  Domain  Measurement  Integrals 

The  estimates  of  the  Karhunen-Loeve  coefficients  omn  are  given  in  terms  of  the  measure¬ 
ments  Z{x,y)  by  equation  (86).  Using  equations  (99)  and  (62)  through  (67)  in  equa¬ 
tion  (86)  yields 


[X~A  fy=B  T 

dmn  =  /  /  Z  ( x,y)Kmn(x,y)dxdy 

J  x—0 


[r=A  fV=B  _  ® 

/  /  Z  (x,  y)  /  .  0mn i0;nni(^i  U )dxdy 

J*= 0  •/»=0  J  =  1 


rx=A  p=B 

!  I  Z  (-X,  y)  {/?mtll  t/'mnl  (x,  y)  +  fimn2lpmnl(x ,  y)  +  0mn3xPmn‘i{x ^  y) 

t=0  J  y=0 


d"/^mn4V,mn4(*^*i  2/)  “f"  ^mn5^mn5(^i  V )  “t*  /^mn60mn6(*^  y)  }  dxdy 
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Tmnl  sinoma;  sin  bny 

0 

0 

Tmn2  sin  amx  sin  bny 

rx=A  ry=B 

=  [  I  ZT(x,y)  < 

Jx= 0  J  y=0 

0mnl 

0 

0 

+  fimn2 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

ft mn3 

Tmn3  sin  amx  sin  bny 

ftmn4 

0 

flmnS 

0 

0 

0 

7mn4  cos  dmx  cos  bny 

0 

0 

Tmn5  sin 

0 

0 

0 

4"  ftmn 6 


0 

0 

0 

0 

0 

7mn6  COS  dmx  sin  bny 


dxdy 


=  /  /  lZi  (xiV)  Z2{x,y)  Z3{x,y)  Z4{x,y)  Z5{x,y )  Z6(x,y)} 

J  X~0  Jy=Q 


X 


ftmnl  Tmnl  sill  Qm£  sin  bny 
ft mn2'ymn2  S1H  dmX  sill  bnXJ 

ft mn3'ymn3  dm  %  sin  bny 
^mn47wn4  COS  dm X  COS  bny 
/?mn57mn5  sin  CLmX  COS  bny 
ftmn6'7mn6  COS  Q mX  sill  bny 


dxdy 


rv~A  rv=B 

=  /  {  ^i(x,y)A»ni Tmnl  sin amx  sin  6ny  +  Z2(x,y)/3mn2ymn2  sin  amx sin  bny 

Jx= 0  J  V=0 
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+Z3(x,  2/)/?mn37mn3  sin  amx  sin  bny  +  Z4(x,  y)^mn4lmM  cos  amx  cos  bny 


+Z$(x,  y)ftmns')mns  sin  am x  cos  bn y  +  Z6(x,y)/3mn6'fmnecosamxsmbny  } 


rx=A  ry=B 

=  fimnilmni  /  /  Zt(x,  y )  sin  amx  sin  bnydxdy 

J  r~Q  J  y=0 


fx=A  ry=B 

+/3mn27mn2  /  /  Z2(x,  y)  sin  amx  sin  bnydxdy 

JrzzO  Jy-0 


rx=A  ry-B 

+/?mn37mn3  /  /  Z3{x ,  y)  sin  amx  sin  bnydxdy 

Jx—O  J  y— 0 


t*-A  ry=B 

+/3mn47mn4  /  /  Z4(x ,  y)  cos  amx  cos  bnydxdy 

J  xzz  0  «/  y=0 


rx—A  ry=B 

+0mn5lmnb  /  /  Zs(x ,  y)  sin  amx  cos  bnydxdy 

Jx=0  *'y=0 


rxzzA  ry=B 

I  I  y)  cos  cimx  sin  6nycfx(iy 

J  x=0  «/y=0 


—  /^mnlTmnl “i“  /^mn2Tmn2^2(^i  ^ ^mn3Tmn3^3(^>  “f"  /^mn4<Tmn4^4(^,i 


"^/^mn5Tmn5^5(^^  ”4“  /^mn6'7mn6^6(^? 


(106) 
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where 


2\(m,n)  = 

J 

rx=A 

1=0  •> 

ry=B 

/  Z\(x,  y)  sin  amx  sin  b„ydxdy 

V=0 

(107) 

rx=A 

ry~B 

11 

s 

Lo  J 

/  Z2(x,  y)  sin  amx  sin  bnydxdy 

y=0 

(108) 

rx=A 

rv=B 

Z3(m,n)  = 

Lo  j 

j  Z3(x,  y)  sin  amx  sin  bnydxdy 

y—O 

(109) 

rx=A 

ry—B 

Z4(m,n)  =  J 

'  I 

x=0  J 

'  Z4(x,  y)  cos  amx  cos  bnydxdy 

y=o 

(110) 

tx—A 

/•v=B 

Z5(m,n)  =  j 

(  I 

r=0  J 

f  Z$(x,  y)  sin  amx  cos  bnydxdy 

v~o 

(1U) 

Z6(m,  n)  -  j 

fx=A 

ry=B 

f  j 
1=0  J 

<  Ze(x,y)  cos  amx  sin  bnydxdy 

y=0 

(112) 

Thus,  it  is  clear  that  if  equations  (107)  through  (112)  can  be  evaluated  and  substituted  in 
equation  (106)  along  with  the  /3mn,’s  from  equation  (105)  and  the  7mm’s  from  equations  (62) 
through  (67)  then  the  estimates  amn  of  the  Karhunen-Loeve  coefficients  can  be  obtained 
from  equation  (106). 

4.9  Signal  Estimates  from  Coefficient  Estimates 

Once  the  Karhunen-Loeve  coefficients  are  estimated,  the  signal  estimates  are  easily  ob¬ 
tained  from  the  basic  representation  of  the  disturbance  potential  T(x,y,  z)  given  by  equa¬ 
tion  (25).  Thus,  at  any  spatial  point  (x,y)  in  the  local  region  defined  by  equation  (1)  and 
any  altitude  h,  the  estimates  of  the  scalar  disturbance  potential,  the  three  components  of 
the  gravity  vector  and  the  six  elements  of  the  gravity  gradients  are  as  given  below: 


Disturbance  Potential  Estimate: 


oo  oo  (  2 

T(x, y,  h)  =  52  52  sin  a™x  sin  ^  ]  /jd' 

m=l  n=l  A  V 


a- 1 


^mnl 


Gravity  Vector  Estimates: 


dT 


(x,y,/i)  =  51  52  “mnc°s<im 
ox  m=1 n=1 


^mn2 


5T 


oo  oo  (  2 

(*.»,»)-  EE  Qmn  sin amx cos  i>ny  | 


I 


&mn3 


dT,  .  x  sr  X2,  -  •  •  i  f-2sgn(/t  + D) 

— (x,  y,h)=  5Z  2J  a"»"  sm  °mX  sin  6"y  | 


Cmn^ 


-Cmn|h+D| 


0mn4 


Gravity  Gradient  Estimates: 


(113) 


xsin&„y  ^=ame  C">"lh+D||  (114) 


6ne-cm„|h+D|l  (115) 
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d2T 


— (x,y,/i)=  52  52®™"smamXsm^ny 


^  ^  2 

V 

^mnS 


(117) 


(^T,  ,  X  ^  ^  -  •  «  -■  f  ~2  U.-Cmn|<.+P|\ 

-^-j(x,  y,M=EE  amn  sm  GmX  sm  6nJ/ 1  y/ABbn  J 


(118) 
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5 


DISCRETE  IMPLEMENTATION 


5.1  Introduction  &  Summary 

In  this  chapter  the  practical  problem  of  discrete  measurements  is  addressed.  The  measure¬ 
ments  are  assumed  to  be  on  a  regular  two-dimensional  grid  in  the  survey  region.  The  data 
samples  are  assumed  to  be  evenly  spaced  in  each  direction  although  there  is  no  restriction 
that  the  sampling  interval  in  either  direction  need  be  the  same.  This  is  particularly  desir¬ 
able  since  the  measurement  samples  in  airborne  gravity  gradiometry  survey  are  obtained 
at  a  much  closer  spacing  (i.e.  1  Km)  on  each  track  even  though  the  spacing  between  the 
tracks  is  further  apart  (i.e.  5  Km).  Since  the  measurement  integrals  involve  sine  and  co¬ 
sine  functions  a  simple  application  of  the  Fourier  transform  formulas  result  in  conversion 
from  integrals  to  summations.  Restricting  the  spatial  frequencies  to  finite  values  these 
double  summations  are  converted  to  simply  left  and  right  transformations  involving  sine 
and  cosine  functions.  Assuming  that  the  number  of  coefficients  in  the  Karhunen-Loeve 
domain  is  equal  to  the  number  of  measurements  in  the  spatial  domain  then  an  obvious 
and  brute  force  approach  to  inverting  the  measurement  integrals  is  by  simply  inverting  the 
sine  and  cosine  transformation  matrices.  The  obvious  drawback  of  this  approach  is  the 
necessity  to  perform  several  large  matrix  inversions.  An  alternate  approach  which  avoids 
this  drawback  altogether  is  presented.  This  approach  is  based  on  the  observation  that  the 
sine  transform  matrices  are  orthogonal  and  that  the  cosine  transform  matrices  also  enjoy 
an  orthogonal  relationship  involving  a  Toeplitz  circulant  matrix  of  alternating  ones  (1) 
and  zeroes  (0).  Utilizing  these  properties  the  continuous  domain  measurement  integrals 
are  easily  discretized  for  the  discrete  grid  of  measurements  without  any  necessity  of  large 
scale  matrix  inversions.  It  should  be  pointed  out  that  by  this  approach  the  transition 
from  the  continuous  domain  to  the  discrete  domain  is  exact  and  not  merely  a  discrete  ap¬ 
proximation  to  a  continuous  algorithm.  After  this,  the  estimation  of  the  Karhunen-Loeve 
coefficients  becomes  a  simple  matter  of  point-by-point  matrix  multiplication  (not  row-by¬ 
column  matrix  multiplication)  and  point-by-point  matrix  summation  of  the  coefficients  of 
the  estimator  gains,  the  discretized  measurement  integrals  and  observation  matrices  asso¬ 
ciated  with  each  of  the  gradient  measurements.  Of  particular  interest,  is  the  fact  that  all 
the  Karhunen-Loeve  coefficients  are  estimated  simultaneously  using  all  the  discrete  two- 
dimensional  grid  measurements  of  all  the  gradients.  Once  the  Karhunen-Loeve  coefficient 
estimates  are  obtained  two-dimensional  grid  estimates  of  the  disturbance  potential,  gravity 
vector  components  and  gravity  gradients  are  easily  obtained  by  point-by-point  multipli¬ 
cation  of  the  Karhunen-Loeve  coefficient  estimates  with  an  observation  matrix  depending 
on  the  type  of  signal  to  be  estimated  and  subsequently  performing  appropriate  left  and 
right  sine  and  cosine  transforms.  As  before,  these  transformations  involve  no  matrix  inver¬ 
sions.  In  addition,  the  two-dimensional  grid  of  signal  estimates  need  not  be  based  on  the 
measurement  grid.  A  finer  interpolated  or  densified  grid  can  be  used  to  obtain  the  grid 
of  signal  estimates.  Also  this  finer  grid  can  be  located  at  any  altitude,  not  necessarily  at 
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the  survey  altitude  by  simply  specifying  the  altitude  parameters  in  the  observation  matrix. 
Thus  downward  continuation  is  performed  automatically. 


5.2  Definition  of  Survey  Grid 

In  order  to  obtain  the  vector  signed  estimates  for  all  the  spatial  points  of  interest  a  definition 
of  the  two-dimensional  grid  is  necessary.  For  the  present  development  the  two-dimensional 
grid  of  data  points  is  defined  to  be  equally  spaced  in  each  direction  with  K  data  points  in 
the  x  -  direction  and  L  data  points  in  the  y  -  direction.  With  Az  and  Ay  the  grid  spacing 
in  the  x  and  y  directions  respectively  the  spatial  domain  D  is  given  as 


(A,  +  l)A*  =  i4  l<k<K 

(123) 

(L  +  l)Ay  =  B  1  <1<L 

5.3  Application  of  Fourier  Transforms 

(124) 

Applying  the  theory  of  Fourier  Transforms  to  the  continuous  domain  measurement  integrals 
given  by  equations  (107)  -  (112)  yields  ( Brigham ,  1973) 

A  OO  OO 

Z\{x,  y)  =  —  /L  ^i(m< n ) sin  “mi  sin  bny 

m= 1  n= 1 

4  00  00 

(125) 

2j(i,  y)  =  n)  sin  amx  sin  bny 

m=l  n=l 

4  00  00 

(126) 

^3 (x,y)  =  —  ]T  Z3(m,n)  sin  amx  sin  bny 

m=l  n= 1 

4  00  00 

(127) 

Z4(x,y)  =  —  £  £  Z4(m,n)  cos  amx  cos  6ny 

m=sl  n=l 

(128) 
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Zs(x,y)  =  ~T5  E  E  Z5(m,n)  sin  amx  cos  bny 

AH  m— 1  n=l 


(129) 


i 


Z6(x,y)  =  -775  £  £  Ze{m,  n)  cos  amx  sin  6ny 

m  =  l  n=l 


(130) 


Limiting  the  spatial  frequencies  to 


1  <m<M 


(131) 


1  <n<N 


(132) 


equations  (125)  through  (130)  become 


4  M  N  _ 

ZtiX'ij)  =  -7-^52  52  Z\(m.n)sin  nmx  sin  6ny 

m=l  n=l 


(133) 


I 


4  M  N  _ 

Z2(x.y)  =  —=  Y,  52  Z2(m,n) sin  a, nxsm  bny 

m=  1  n=I 


(134) 


4  M  N  _ 

Z3(x,y)  =  75^1  Z3(rn,  n )  sin  rtmxsm6n  y 

m  —  1  Ti— 1 


(135) 


4  «  w  . 

Z4(x,y)  =  ~Tn  52  52  z4(m,n)  cos  nmx  cos  bny 
AH  m=1  n=1 


(136) 


4  At  N 

Z$(x,y)  =  —  EE  Z5(m,  ?) )  sin  «mx  cos  5„y 

m=l  n=  1 


(137) 
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(138) 


4  M  N 

y)  =  jg  Y  Y  Zeim,  n)  cos  amx  sin  bny 

m=l  n=  1 


Discretizing  as  per  equations  (123)  and  (124)  such  that 


x  =  xk  \<k<K  (139) 


y  =  yi  1  <1<L  (140) 


and  writing  in  matrix  form  equations  (133)  through  (138)  appear  as 


[Z,]  =  _4_  [SAX]  [z,]  [SAF]r 
KxL  AB  KxMMxNNxL 


(141) 


[z2]  =  A  i5A*i  [5AFir 

KxL  AB  KxMMxNNxL 


(142) 


[Z3]  =  [SAX]  [z3]  [SAF]r 
KxL  AB  KxMMxNNxL 


(143) 


[Z4]  =  [C AX]  [z4]  [CAF]r 

KxL  AB  KxMMxNNxL 


(144) 


[Z5]  =  _4_  [SAX]  [Z5J  [CAY}t 
KxL  AB  KxMMxNNxL 


(145) 
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(146) 


[Ze]  =  J_  [C AX]  [Z6]  [5AF]7 
KxL  AB  KxMMxNNxL 


where 

K,L  -  Number  of  measurements  in  spatial  domain  x,y  plane 
M.  N  Number  of  coefficients  in  Karhunen-Loeve  domain 


Zi(xi , y 1 ) 

Zi(xi,y2)  ■ 

Z\(xuVl) 

[2,]  _ 
KxL 

Zi(.r2,yi) 

.  Zx(xK,y ,) 

■  Z\{xx,yi)  . 

Z2(xx ,  j/i ) 

Z2(xi,y2)  ■ 

■  Z2{xx,yi) 

[^]  _ 
KxL 

Z2(x2,yx ) 

.  Z2(xK,yi) 

■  Z2(xK,yi)  . 

Z3{xi,yi) 

Z3(xuy2)  ■ 

■  Z3(xx,yi) 

[Z3)  _ 
KxL 

Z3{x2,yi) 

.  Z3(xK,t/i) 

■  Z3(xK,yc)  . 

Za(x\ ,  j/i ) 

Z4(xily2)  ■ 

•  Z4{xx,yi) 

\z*]  _ 

KxL 

Z4(x 2,  y\ ) 

.  Z4(xh-,yt ) 

... 

■  Z4{xK,yL)  _ 

(147) 


(148) 


(149) 


(150) 
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r 


f 

i 

i 

> 


f 

k 


► 


I 


[Zs]  _ 
K  x  L 


[2«]  = 

K  x  L 


Zs(xi,y\) 

Z&(x2,yi) 

Zb(x\,y2)  ■ 

•  Zs(xuyL) 

(151) 

.  Zs(xK,yi) 

... 

•  Zs{xK,yL)  . 

'  Z6(xi,yj) 

Z6(x2,yi) 

Ze(xi,  y2)  ■ 

■  Ze{xi,yL) 

(152) 

.  Ze(xK,yi) 

... 

•  Ze(xx,yL )  . 

[SAX] 

KxM 


[CAX] 

KxM 


[say]t 

N  x  L 


[CAYjT 
N  x  L 


sin  7r(  1  x  1)^2C 
sin  7r(2  x  1)^ 

sinir(lx2)^f  •• 
sin  ir(2  x  2)^ 

•  sin  ir(  1  x  Af)M.  ‘ 

•  sin  ir(2  x  M)~ 

(153) 

.  sin7r(A'  x  1)^ 

sinir(A  x  2)~  •  • 

■  sin  7r( A  x  M)^f  . 

COS7r(l  X  1)^- 
COS  7r(  2  X  1)A^ 

COS7r(l  x  2)^  • 

cos  ir(2  x  2)^j-  • 

•  cos7r(l  x  M)^r 
■  COS7t(2  x  A/)  — 

(154) 

COS  7r(/\  X  1)^ 

cos  ir(A'  x  2)—^-  • 

■  cos  it  (A"  x  M)^~ 

sin7r(l  x  1)^- 
sin  ?r( 2  x  1)^- 

sinir(l  x  2)^  • 

sin  7r(2  X  2)^- 

sin  ?r(  1  x  A)^ 

•  •  sin  ir(2  x  L)^g- 

(155) 

sin  7r(]V  x  1)^- 

sin  n(N  x  2)^g-  • 

•  •  sin  ir(  N  x  A )^f- 

cosir(l  x  1)A£ 
COS7r(2  x  1)^ 

cos  7r(  1  x  2)^ 
cos7r(2  x  2)^  • 

cos  7r(l  x  A)^ 
cos?r(2  x  L)^f- 

(156) 

.  cosir (N  x  1)^- 

cos  7t(AT  x  2)^p  • 

■  ■  COS  7T (N  x  A)^ 
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Thus, 


[SAX(i,j)]  =  smn(ixj)^£-  \<j<M  (157) 

[SAF(i,  j)]T  =  sin  ir(i  x  j)*£  1  <  ;'  <  (158) 

[CAX(iJ)\  =  cost (i  x  |  *  (159) 

[CAF(i,»]T  =  cos  7r( ?  x  j)^  }  |  *  (160) 


5.4  Brute  Force  Approach  Requiring  Matrix  Inversions 

Assuming  that  the  number  of  coefficients  in  the  Karhunen  Loeve  domain  is  equal  to  the 
number  of  measurements  in  the  spatial  domain  satisfying  the  conditions 

M  =  K  (161) 


N  =  L  (162) 

then  an  obvious  and  brute  force  approach  to  inverting  equations  (141)  through  (146)  is  by 
simply  inverting  the  matrices  such  that 


Z1]  =  ^(5AXr,[Z1]((5Ar]T)'1 


(163) 
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Z2]  =  ~  [SAX]'1  [Z2]  ([5AF]t)_1 


(164) 


lZ3]  =  ^  [SAX]"1  [Z3]  ([SAy]7’)'1  (165) 


[Z4j  =  [C AX]'1  [Z4]  ([CAr]r)_1  (166) 


[is]  =  ~  [SAX]-1  [Zs]  ([CAy]7")"1  (167 ) 


[i«l  =  [CAX]-1  [Ze]  ([SAy]7-)-1  (168) 

The  obvious  drawback  of  this  approach  is  the  necessity  to  perform  several  large  matrix 
inversions-  In  the  next  section  an  alternate  approach  is  outlined  which  avoids  this  drawback 
altogether. 

5.5  Alternate  Approach  Avoiding  Matrix  Inversions 

Bose  et.  al.  (1988)  have  recently  shown  that  the  particular  matrices  under  consideration 
here  given  by  equations  (153)  through  (156)  enjoy  the  following  properties  outlined  below: 
1)  Sine  Transform  matrices  given  by  equations  (153)  and  (155)  satisfy  the  condition 


[SAX]r  [SAX]  =  K  +  1  ,f 

K  x  K  K  x  K  2  ‘  1 


(169) 


where  [I]  is  a  K  x  K  identity  matrix.  The  above  condition  leads  to 


([SAXHSAX])-1  =  ^—[7] 
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(170) 


[SA*]-'  = 

2)  Cosine  transform  matrices  given  by  equations  (154)  and  (156)  satisfy  the  condition 


[cax]t[caa:]  _  k  +  i 

KxK  KxK  ~  2 


[I\  ~  [C] 


(171) 


provided  K  is  even  and  not  odd  and  [C]  is  a  K  X  K  Toeplitz  Circulant  matrix  of  alternating 
one’s  (1)  and  zeroes  (0),  a  4  x  4  example  of  which  is 


[Cj  = 


10  10 
0  10  1 
10  10 
0  10  1 


(172) 


The  properties  of  Toeplitz  Circulant  matrices  are  such  that 


(^~ [/]  -  [C])"1  =  ~ ^([/]  +  2[C])  K  -  even 


(173) 


which  leads  to 


([CAXHCAX])-1  =  +  2  [C]) 


or 


fCAX]"1  =  ^-t([/]  +  2[C])[CA.Y]-1  (174) 

Using  the  assumptions  given  by  equations  (161)  and  (162)  and  the  properties  given  by 
equations  (170)  and  (174)  and  the  domain  definition  given  by  equations  (123)  and  (124) 


in  equations  (141)  through  (146)  yields 


[Zil  =  AB  [SAX]r[Z1][5Ay^] 

K  x  L  (K  +  1  )(L  +  1)  KxK  KxL  LxL 


(175) 


[Z2]  =  AB  [SAX]r[Z,][SAr] 

KxL  (A+1)(I  +  1)  KxK  KxL  LxL 


[Z3\  =  AB  [5AJT]r(Z3J[5Ar] 

KxL  (K  +  1)(L  -f  1)  KxK  KxL  LxL 


[Z4]  =  AB  [/  +  2C][CAX]t[Z4][CA>'][7  +  2  C] 

KxL  (A' +!)(£  + 1)  KxK  KxK  KxL  LxL  LxL 


[Z5]  =  AB  [SAJr)7’(Z5][CAF][/-f  2C] 

KxL  (A'  +  1){L  +  1)  KxK  KxL  LxL  LxL 


[Z6]  =  AB  [I  +  2C][CAX]r[Z6][5AK] 

KxL  (A'+1)(A  +  1)  KxK  KxK  KxL  LxL 


(ISO) 


5.6  Discrete  Estimation  of  Karhunen-Loeve  Coefficients 

Representing  ®  as  point- by- point  matrix  multiplication  and  not  row-by-column  matrix 
multiplication,  equation  (106),  under  the  assumptions  given  by  equations  (161)  and  (162) 
can  be  written  in  discrete  matrix  form  as 


[«mu]  _  /?mrn  ®  ^i(m,n)  @7mnl  -f  0mni  ®  Z2(m,n)  ®  ymn2 
KxL  KxL  KxL  KxL  KxL  KxL  KxL 


56 


(181) 


+  0mn 3  ®  Z3(m ,  n)  (8)  7m n3  +  0mn4  ®  Z4 (m,  Tl)  <8>  7m n4 
KxL  KxL  KxL  KxL  KxL  KxL 


+/?mn5  ®  Zs(m,  n )  ®  7mn5  +  ^mn6  ®  ^6(m,  n)  ®  Jmn6 
KxL  KxL  KxL  KxL  KxL  KxL 


where  the  0mni  s  are  given  in  equation  (105),  7m„i’s  are  given  in  equations  (62)  through  (67) 
and  the  Zi(m,n)' s  are  given  in  equations  (175)  through  (180). 

5.7  Two-dimensional  grid  estimates  of  signals 

From  equations  (113)  through  (122)  it  is  obvious  that  the  signal  estimates  are  available  at 
any  altitude  h.  However,  since  the  estimates  of  the  Karhunen-Loeve  coefficients  are  given 
by  equation  (181)  wherein  amn  is  available  in  a  K  x  L  matrix  form,  therefore,  in  order  to 
obtain  the  signal  estimates  the  conditions  given  by  equations  (161)  and  (162)  are  necessary. 
But  the  two-dimensional  grid  of  signal  estimates  need  not  be  based  on  the  measurement 
grid  defined  by  equations  (123)  and  (124).  A  finer  interpolated  grid  defined  by 


(P +  !)£*  =  A  SAT  <  A.Y  (182) 


(Q  +  l)Y,Y  =  B  HY  <  AY  (183) 


can  be  used  to  obtain  the  grid  of  signal  estimates.  Thus,  utilizing  this  grid  definition 
along  with  necessary  conditions  given  by  equations  (161)  and  (162)  in  equations  (113) 
through  (122)  yields  the  following  two-dimensional  grid  of  signal  estimates. 

Disturbance  Potential  Estimate: 


[T]  =  [5E.Y]([amn]  ®  [fUi])[SS5T 

PxQ  PxKKxL  KxLLxQ 


(184) 


Gravity  Vector  Estimates: 
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[f  ]  =  [<?£jr]([smn]  ®  [^mn2])[5£y]] 

PxQ  PxKKxL  KxLLxQ 


f  ]  =  [SEX]([am„]  0  [0mn3])[C£y]7 

PxQ  PxKKxL  KxLLxQ 

f]  =  [S£X]([amn]  0  [0mn4])[S£y]7 

PxQ  PxKKxL  KxLLxQ 
Gravity  Gradient  Estimates: 


[fg]  =  [5£A'j((amnj  0  [0m„5J)[$£yia 


PxQ 

P  x  K  K  x  L 

K  x  L  L  x  Q 

an  _ 

dyi 

[S£jr]([smn]  0  [0mn6])[ssy]r 

(ISO) 

PxQ 

P  X  K  K  X  L 

K  x  L  L  x  Q 

an  _ 

dz^ 

[5EX]([am„] 

®  [0mn7])[5Sy]T 

(190) 

PxQ 

P  x  K  K  x  L 

K  x  L  L  x  Q 

an  _ 

dxdy 

[CEX]([aTOn 

®  [^n8])[csy]r 

(191) 

PxQ 

P  x  K  K  x  L 

K  x  L  L  x  Q 
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(192) 


f 


d2T 

dydz 


=  [S£X]([Smn]  ®  [0mn9])[CEy]T 


PxQ  PxKKxL  KxLLxQ 


d2T 


PxQ 


[C£X]([5mn]  ®  [0mnlo])[SEF]T 

PxK  KxL  KxL  LxQ 


(193) 


where  ®  has  been  defined  earlier  as  point-by-point  matrix  multiplication,  0mnt's  are  defined 
in  equations  (113)  through  (122)  and 


[SEX] 

P  x  K 

sin7r(l  x  l)^p 
sin  7r(2  x  1)~X 

sin7r(l  x  2)^-  • 

sin7r(2x2)^f  ■■ 

•  sin7r(l  x  K)^~- 
■  sin  tt(2  x  K)^~- 

(194) 

.  sin  7r(P  : 

xi)2£ 

sin  7r(P  x  2)^  ■  ■ 

•  sin7r(P  x  A')Sf  . 

[CSX] 

PxK 

cos7r(l  x  l)^r 
cos  7r(2  X  1)-^ 

COS  7f(  1  x  2)^ 
cos7r(2  x  2)~  ■ 

■  •  COS7r(l  X  K ) “ 

••  cos7r(2xA")^ 

(195) 

- 

cos  7r(P  ) 

<D¥ 

cos  t r(P  x  2)^  • 

•  •  cos  7 r(P  x  K)^j- 

[SEy]r 

LxQ 

sin7r(l 

sin7r(2 

X  X 

sin  7r(  1  x  2)^  ■ 
sin  7r(2  x  2)^  • 

■  •  sin  7r(l  x  Q)^f 

■  ■  sin  7r(2  x  Q)^§- 

(196) 

.  sin  tt(L 

xi)¥ 

sin  n(L  x  2 )%g-  ■ 

■  ■  sin  7 t{L  x  Q)^j-  . 

[csyf 

LxQ 

COS7t(1 
COS  7r(2 

xl)f 

xl)# 

sin7r(l  x  2)^  • 

sin7r(2  x  • 

■  ■  cos7r(l  x  Q)^-s- 

■  •  sin7r(2  x  Q )~ 

(197) 

.  cost r(L 

xl)f 

sin  7 r(L  x  2)^-  ■ 

■  •  sin  7 r(L  x  Q)^jf 
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6  CONCLUSIONS 


6.1  Summary  of  Research  Performed 

A  need  to  devise  a  methodology  to  process  two-dimensional  grids  of  gravity  gradients  at 
survey  altitude  to  yield  gravity  disturbance  vector  estimates  at  the  surface  of  the  earth, 
motivated  this  research.  The  measured  gravity  gradients  are  the  six  elements  of  the  gra¬ 
dient  tensor.  An  actual  airborne  gradiometer  survey  over  an  area  300x300  Km  consisting 
of  bidirectional  flight  paths  5  Km  apart  with  along-track  sampling  intervals  of  1  Km  will 
result  in  approximately  220,000  measurements.  The  main  problem  with  the  determination 
of  the  gravity  field  from  airborne  gradiometry  is  the  huge  amount  of  gradient  data  collected 
during  a  gradiometry  survey.  The  primary  objective  of  this  research  effort  was  to  solve 
the  problem  of  processing  all  the  airborne  gravity  gradient  measurements  simultaneously 
in  a  computationally  efficient  manner  without  neglecting  gradiometer  measurement  noise. 

The  approach  taken  here  is  to  exploit  the  marriage  of  physical  theory  of  geodesy  and 
random  process  theory.  The  gravity  signal  model  used  is  obtained  by  solving  Laplace’s 
equation  with  the  unknown  mass  distribution  below  the  surface  of  the  earth  modelled  as 
a  two-dimensional  white  noise  layer  representing  the  vertical  derivative  of  the  disturbance 
potential  to  any  pre-specified  order.  This  results  in  a  series  solution  of  the  disturbance 
potential  wherein  the  unknown  coefficients  of  the  expansion  are  forced  to  be  uncorrelated 
by  invoking  the  Karhunen-Loeve  condition.  This  resulting  disturbance  potential  covariance 
obtained  from  this  model  is  both  non-stationary  and  non-isotropic. 

The  six  (6)  gravity  gradients  being  functionals  of  the  disturbance  potential  were  rep¬ 
resented  in  terms  of  the  Karhunen-Loeve  series  expansion  of  the  disturbance  potential 
resulting  in  six  basis  functions.  These  basis  functions  were  shown  to  be  orthogonal.  The 
measurement  model  was  chosen  to  be  these  six  gravity  gradients  in  the  survey  region  con¬ 
taminated  by  additive  white  noise.  The  estimation  problem  was  to  take  all  the  six  gravity 
gradient  measurements  in  the  survey  region  and  obtain  estimates  of  the  gravity  vector 
components  at  the  surface.  This  estimation  problem  was  shown  to  be  equivalent  to  simply 
estimating  the  Karhunen-Loeve  coefficients  from  all  the  gradient  measurements. 

Under  the  assumption  of  gaussianness  of  noise  statistics  the  optimal  estimator  was 
represented  by  expressing  each  Karhunen-Loeve  coefficient  as  a  linear  functional  of  all 
the  measurement  data.  Each  of  these  linear  functionals  was  specified  by  an  associated 
weighting  function.  The  form  of  each  weighting  function  was  found  by  employing  the 
orthogonality  principle  of  linear  mean  square  estimation  theory. 

The  practicality  of  discrete  measurements  motivated  discretization  of  the  continuous  al¬ 
gorithm.  Only  equally  spaced  discrete  points  were  considered.  Continuous  two-dimensional 
measurement  integrals  were  converted  to  summations  by  utilization  of  special  orthogonal¬ 
ity  properties  of  sine  and  cosine  transforms  discovered  during  the  course  of  this  research. 
The  discrete  algorithm  turned  out  to  be  expressible  as  a  sequence  of  matrix  multiplications 
without  any  necessity  of  matrix  inversions  whatsoever.  It  was  shown  that  for  evenly  spaced 
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data  points  the  discrete  algorithm  derived  is  exact  and  not  merely  an  approximation  to  a 
continuous  algorithm. 

Once  the  Karhunen-Loeve  coefficient  estimates  were  obtained  from  discrete  data,  two- 
dimensional  grid  estimates  of  the  disturbance  potential,  gravity  vector  components  and 
gravity  gradients  were  easily  obtained  by  point-by-point  multiplication  of  the  Karhunen- 
Loeve  coefficient  estimates  with  an  observation  matrix  depending  on  the  type  of  signal 
to  be  estimated  and  subsequently  performing  appropriate  left  and  right  sine  and  cosine 
transforms.  As  before,  these  transformations  involve  no  matrix  inversions.  In  addition, 
the  two-dimensional  grid  of  signal  estimates  need  not  be  based  on  the  measurement  grid. 
A  finer  interpolated  or  densified  grid  can  be  used  to  obtain  the  grid  of  signal  estimates. 
Also  this  finer  grid  can  be  located  at  any  altitude,  not  necessarily  at  the  survey  altitude 
by  simply  specifying  the  altitude  parameters  in  the  observation  matrix.  Thus  downward 
continuation  is  performed  automatically. 

6.2  Highlights  of  Research  Achievements 

A  methodology  of  post-mission  data  processing  for  the  gravity'  gradiometer  survey  system 
is  presented.  The  highlights  of  the  technique  include: 

1.  The  model  is  derived  from  the  physical  theory  of  geodesy  and  is  not  based  upon 
empirical  assumptions  of  correlation  functions  or  power  spectral  densities. 

2.  Model  can  accommodate  multiple  two-dimensional  white  noise  layers  below  the  sur¬ 
face  of  the  earth. 

3.  Each  layer  can  model  the  vertical  derivative  of  the  disturbance  potential  to  any  order. 

4.  Non-zero  boundary'  values  for  the  disturbance  potential  on  the  exterior  of  the  survey' 
region  permitted. 

5.  The  model  is  such  that  at  any  given  spatial  point  the  gravity  field’s  correlation  with 
neighboring  points  is  preserved  and  correlation  in  any  direction  is  not  ignored. 

6.  The  correlation  of  the  gravity  field  with  increasing  distance  is  not  ignored. 

7.  The  estimation  algorithm  does  not  enforce  any  unnecessary  limitation  of  causality 
on  the  data  inasmuch  as  no  one-dimensional  scanning  is  performed. 

8.  Each  coefficient  in  the  series  expansion  for  the  gravity  field  is  estimated  using  all  the 
six  inline  and  crossline  gravity  gradiometer  data  simultaneously. 

9.  The  estimation  algorithm  can  handle  gradient  data  given  in  two-dimensional  grids 
at  the  same  or  different  altitudes  on  or  above  the  surface  of  the  earth. 


10.  Downward  continuation  of  the  gravity  field  from  measurements  at  the  survey  height 
above  the  surface  of  the  earth  is  automatically  done  without  any  loss  of  accuracy. 

11.  Interpolation  of  estimates  between  grid  measurements  performed  automatically  by 
employing  a  denser  grid  for  prediction. 

12.  Different  apriori  accuracies  can  be  assigned  to  measurements  from  different  gradiome- 
ter  inline  and  crossline  measurements. 

13.  Correlated  noise  sources  can  be  accommodated  to  the  extent  that  they  cam  be  rep¬ 
resented  by  the  basis  functions. 

14.  The  estimation  algorithm  requires  no  matrix  inversions. 

15.  The  discrete  algorithm  is  exact  and  not  merely  an  approximation  to  a  continuous 
integral. 

16.  Measurement  data  must  be  in  planar  gridded  form. 

6.3  Recommendations  for  Future  Research 

Listed  below  are  suggestions  for  future  work  to  enhance  this  method  of  post-mission  data 
processing  for  the  gravity  gradiometer  survey  system: 

1.  The  White  Noise  Layer  (WNL)  model  used  in  the  algorithm  is  based  on  a  single  layer 
of  white  noise  below  the  surface  of  the  earth.  Multiple  layers  can  be  accommodated 
and  should  be  investigated  to  examine  its  effect  on  the  modeling  sensitivity. 

2.  Other  sensors  particularly  those  providing  long  wavelength  information  should  be 
incorporated  in  the  estimation  algorithm. 

3.  Measurement  error  models  more  elaborate  than  the  present  simple  white  noise  model 
may  be  worth  investigating  to  better  model  the  measurement  errors. 

4.  Develop  implementable  algorithms  to  take  care  of  the  non- zero  boundary  conditions 
in  a  manner  consistent  with  survey  data. 

5.  Examine  the  relationship  of  this  estimator  with  that  of  Wiener  filtering,  Fourier 
transforms  and  Least  Squares  Collocation. 

6.  Develop  algorithms  to  obtain  a  theoretical  error  covariance  by  using  Kronecker  matrix 
products  to  the  discrete  algorithm. 

7.  Use  error  analysis  algorithms  to  perform  combination  solution  trade-off  analysis. 
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8.  Investigate  the  performance  of  the  algorithm  for  modeling  error  both  for  signal  model 
and  for  noise  models. 

9.  Investigate  optimization  of  the  discrete  algorithm  from  the  standpoint  of  minimizing 
number  of  operations  as  is  done  in  Fast  Fourier  Transform  techniques. 

10.  Exploit  parallel  processing  techniques  to  implement  the  algorithm  by  means  of  special 
purpose  microprocessors  specifically  designed  for  this  algorithm. 

11.  Investigate  applications  of  the  algorithm  to  fields  other  than  geodesy  such  as  mag¬ 
netics,  atmospheric  sciences,  topography,  image  processing,  etc. 
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